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ABSTRACT 


Exact,  efficient  solutions  for  realistic  queueing  network 
models  of  computer  systems  are  hard  to  achieve.  Many 
approximate  solution  techniques  have  appeared  in  the  last 
few  years.  One  of  the  most  important  problems  in  using  an 
approximate  solution  technique  is  the  error  behaviour  of 
that  approximation.  The  lack  of  proper  error  estimates  for 
most  approximations  makes  users  of  these  approximations 
rightfully  skeptical.  This  thesis  addresses  approximate 
solution  techniques  and  their  error  behaviours.  Three  major 
contributions  are: 


First,  different  ways  of  improving  an  approximation  are 
examined.  A  unified  view  for  various  approximations  is 
presented.  The  error  introduced  in  an  approximation  is 
divided  into  two  classes:  global  and  local.  Based  upon  this 
classification  of  errors,  examples  are  given  to  show  that 
the  unified  view  can  be  used  to  improve  some  of  the  existing 
approximations. 

Second,  the  unified  approach  also  helps  in  defining  good 
approximations.  We  propose  a  new  approximate  solution 
technique  called  generalized  product  form  (GPF).  Many  of 
the  existing  solution  techniques  can  be  viewed  in  the 
framework  of  GPF.  Two  realizations  of  the  GPF  are  discussed 
and  their  accuracies  are  examined  with  the  help  of  examples. 
The  GPF  framework  should  be  useful  in  introducing  new 


approximations  and  in  analyzing  the  accuracy  of  an  existing 
approximation. 

Third,  experiments  examining  the  error  behaviour  of  various 
approximations  are  presented.  The  configurations  of  the 
queueing  networks  for  the  experiments  are  chosen  based  on 
some  actual  models  of  computer  systems.  Results  obtained  in 
these  experiments  are  analyzed  and  this  analysis  leads  to 
new  insights  into  the  accuracy  of  the  approximations. 
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Chapter  1 


INTRODUCTION 


Impressive 
to  be  made,  with 
increasing.  These 
by  new  applications 
expectations.  An 
heavily  dependent 
systems . 


advances  in  computing  hardware  continue 
hardware  cost  decreasing  and  speed 
advances  have  been  outstripped,  however, 
of  computers  and  ever  increasing  user 
increasing  number  of  organizations  are 
on  efficient  and  reliable  computer 


The  large  investment  in  computing  systems  as  well 
as  their  size,  complexity,  and  changing  nature  all 
contribute  to  the  importance  of  computer  performance 
evaluation.  The  demand  for  computing  power  has  to  be 
forecast,  and  future  resource  requirements  have  to  be 
estimated.  Performance  evaluation  is  an  ongoing  task. 
Tuning  the  system  to  improve  its  performance  and  predicting 
future  needs  and  performance  are  just  two  examples  of 
performance  evaluation. 

It  is  impossible  to  analyze  a  computer  system  in 
its  full  complexity.  Instead  the  modeling  process  allows 
us  to  concentrate  on  only  those  aspects  of  the  system  that 
are  relevant  to  our  particular  concerns.  Models  are  easier 
to  analyze  since  they  are  simpler  in  structure  than  the 
systems  they  represent. 
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A  computing  system 

is 

composed 

of 

multiple 

resources  such  as  the  CPU, 

I/O 

channels , 

disks  , 

and  main 

memory.  Most  large-scale  computer  systems  are 
multiprogramming  systems;  that  is,  different  programs 
(customers)  can  receive  service  simultaneously  at  different 
computer  system  resources.  For  example,  while  one  program 
is  executed  at  the  CPU,  another  one  can  receive  service  at 
one  of  the  channels.  It  is  natural  to  represent  a  computing 
system  in  a  model  as  a  network  of  resources  (service 
centres)  where  jobs  can  move  probabilistically  about  the 
network  from  service  centre  to  service  centre.  Since  more 
than  one  customer  may  require  service  at  a  particular 
service  centre  at  the  same  time,  queues  of  waiting  customers 
naturally  arise. 

Network  of  queues  models  are  generally  solved  by 
simulation  or  analysis.  Simulation  has  been  a  standard 
performance  evaluation  tool  and  is  the  most  general,  most 
flexible,  and  most  powerful  technique  for  studying  and 
predicting  system  performance.  Due  to  high  cost  and 
uncertain  statistical  accuracy,  however,  simulation  is  not 
always  the  first  choice  of  systems  analysts. 

Analytic  models  of  computer  systems  have  proven  to 
be  cost  effective  tools  for  evaluating  computer  systems 
[Graham  1978].  For  example,  they  have  been  used 
successfully  to  evaluate  alternative  proposals  for  computer 
configurations  [Bhandiwad  and  Williams  1975,  Sevcik  and 
Unruh  1975],  to  evaluate  alternative  design  decisions 
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[Lassettre  and  Scherr  1972,  Browne  et  al.  1975],  and  to 
evaluate  the  performance  of  an  existing  system  [Lipsky  and 
Church  1977,  Buzen  1978].  The  use  of  queueing  network 
models,  however,  depends  upon  the  availability  of  efficient 
solution  algorithms  to  calculate  the  performance  measures  of 
interest. 

Queueing  network  models  can  be  classified  into 
three  categories,  according  to  their  solvability  [Chandy  and 
Sauer  1978]:  tractable  models  -  those  that  can  be  solved 
exactly  and  efficiently;  intolerably  slow  -  those  models 
that  can  be  solved  exactly,  but  cannot  be  solved 
efficiently;  and  intractable  models  -  those  that  do  not  have 
known  exact  solutions.  In  the  remainder  of  the  thesis, 
unless  otherwise  stated,  we  consider  only  efficient  solution 
algorithms  to  intolerably  slow  models.  The  solutions 
obtained  by  using  an  efficient  solution  algorithm  may  be 
exact  or  approximate. 

Exact,  efficient  solution  algorithms  exist  only  for 
simple  network  models.  As  the  amount  of  detail  in  a 
queueing  network  model  increases  to  reflect  more  and  more 
features  of  the  real  system,  the  complexity  of  the  model 
drastically  increases.  Not  only  might  an  exact  solution 
not  exist,  but  even  an  approximate  solution  might  prove  to 
be  intolerably  slow.  Thus,  there  is  a  tradeoff  involved 
between  obtaining  an  exact  solution  to  a  simple  model  and 
obtaining  an  approximate  solution  to  a  more  realistic  model. 
According  to  Kleinrock  [1976]: 
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" It  is  not  hard  to  convince  oneself  that  queueing  theory  is 
rather  difficult  and  that  exact  results  are  hard  to  obtain; 
in  fact  many  of  the  interesting  phenomena  have  not  yet 
yielded  to  exact  analysis  (and  perhaps  never  will)." 

Since  exact  solutions  to  more  complicated  networks 
are  hard  to  achieve,  we  must  improve  existing  approximate 
solution  algorithms  or  introduce  better  approximation 
techniques.  Many  approximation  techniques  have  appeared  in 

the  last  few  years  based  on  diverse  ideas  from  different 
fields,  such  as  Norton's  theorem  in  electrical  networks 
[Chandy  et  al.  1975a]  and  near -decomposabil i ty  in  economic 
systems  [Courtois  1977],  In  Chapter  2,  we  shall  discuss 
these  and  other  approximate  solution  techniques  in  some 
detail . 

One  of  the  most  important  problems  in  using  an 
approximate  solution  technique  is  its  error  behaviour.  The 
lack  of  proper  error  estimates  for  most  approximations  makes 
users  of  these  approximations  rightfully  skeptical.  This  is 
especially  a  problem  because  of  the  diverse  origins  of  the 
different  approximations.  Not  only  are  the  error  behaviours 
of  the  different  approximations  quite  different,  but  their 
error  behaviours  have  only  been  studied  individually  (if  at 
all)  . 

This  thesis  addresses  approximate  solution 
techniques  and  their  error  behaviour, 
principal  contributions  of  this  thesis: 


There 


are  three 
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(1)  A  unified  approach  to  understanding  and  improving 
different  approximate  solution  techniques  is  presented. 

(2)  An  approximate  solution  technique  that  shows  better 

accuracy  than  most  of  the  existing  ones  is  proposed. 

(3)  The  error  behaviour  of  existing  approximate  solution 
techniques  is  analyzed  using  experiments.  These  experiments 
suggest  further  ways  to  improve  existing  approximate 
solution  techniques. 

The  next  three  sections  outline  the  importance  of  each  of 
these  contributions  in  computer  system  modeling.  The  final 
section  of  this  chapter  describes  the  organization  of  the 
remainder  of  the  thesis. 

1.1  Unified  Approach  to  Improving  Approximations 

Approximations  are,  in  general,  intuitive 
heuristics,  and  room  for  improvements  still  exists.  In  the 
past,  several  attempts  have  been  made  to  improve  individual 
approximations.  For  example,  Sauer  and  Chandy  [1975] 

improved  the  Norton's  theorem  reduction  proposed  in  [Chandy 
t  al.  1975a].  Sevcik  et  al.  [1977]  analyzed  Norton's 
theorem  reduction,  as  modified  by  Sauer  and  Chandy,  in  more 
detail  and  suggested  improvements.  Gelenbe  and  Pujolle 

[1976]  and  Yu  [1977]  attempted  to  improve  the  performance  of 
diffusion  approximations  suggested  by  Kobayashi  [1974]. 

We  shall  show  that  most  approximate  solution 
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techniques  have  common  sources  of  errors  and  that  there  is  a 
unified  approach  to  improving  many  of  them.  Existing 

approximations  have  a  common  simplification  approach 
decomposition.  The  sources  of  error  in  different 

approximations  such  as  Norton's  theorem  reduction,  extended 
product  form,  and  diffusion  can  be  categorized  as  being 
either  local  or  global.  We  shall  show  that  some  of  the 
recent  approaches  to  improving  Norton's  theorem  reduction 
can  also  be  used  to  improve  diffusion  approximations. 


1.2  A  New  Approximation  Technique 


Considerable  effort  is  being  devoted  to  the  problem 
of  developing  new  approximate  solution  techniques.  For 
example,  Shum  and  Buzen  [1977]  proposed  an  extended  product 
form  and  Marie  [1979]  suggested  an  iterative  approximate 
solution  technique.  Improvements  to  the  existing 
approximations,  however,  can  be  made,  and  more  accurate 
approximate  solution  techniques  are  still  needed. 

The  unified  approach  presented  here  also  helps  in 
defining  better  approximations  than  most  existing  ones.  We 
shall  propose  a  new  approximate  solution  technique  called 
generalized  product  form.  Many  of  the  existing 
approximations  can  be  viewed  in  the  framework  of  this  new 
approximation.  We  shall  show  that  approximations  based  on 


this  new  generalized  product  form  technique  have  better 
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accuracy  than  most  existing  approximations. 


1.3  Error  Behaviour 


A  user's  confidence  in  an  approximate  solution 
technique  can  only  be  based  upon  the  error  behaviour  of  the 
technique.  In  only  a  few  cases,  however,  is  any  formal 
error  analysis  available  for  the  approximations  (e.g., 
[Courtois  1977,  Yu  1977]).  Analytic  error  bounds  are  very 
hard  to  obtain  [Chandy  and  Sauer  1978],  and  we  are  forced  to 
rely  upon  empirical  studies  and  intuition. 

This  thesis  presents  experiments  examining  the 
error  behaviour  of  various  approximate  solution  techniques. 
The  parameters  and  the  configurations  of  the  queueing 
networks  chosen  for  the  experiments  are  based  on  actual 
models  of  computer  systems  as  well  as  theoretical  results 
developed  in  this  thesis.  These  experiments  give  error 

behaviours  of  various  approximations,  which  provides  insight 
into  the  problems  of  using  these  approximations. 

1.4  Organization  Of  The  Thesis 

Chapter  2  presents  an  overview  of  efficient 
solution  techniques  of  aueueing  network  models.  It 
discusses  sufficient  conditions  for  a  network  to  have  an 


exact  solution,  the  implications  of  these  conditions  on  the 
modeling  power  of  queueing  networks,  the  need  for 
approximate  solution  techniques,  and  various  approximate 
solution  techniques. 

Chapter  3  presents  a  unified  approach  to  improving 
approximations.  Based  on  classifying  the  various  sources 
of  errors  into  global  and  local  errors,  this  chapter  unifies 
different  approximate  techniques  in  terms  of  the  sources  of 
errors  present  in  them  and  shows  that  these  unifying 
criteria  can  be  further  used  to  improve  the  accuracy  of  some 
approximation  techniques. 

In  Chapter  4,  we  present  an  approximate  solution 
technique  called  generalized  product  form.  The  formulation 
of  this  approximation  is  based  upon  the  unified  approach 
presented  in  Chapter  3.  Two  realizations  of  generalized 
product  form  (GPF) ,  GPF  with  diffusion  and  GPF  with 
simulation,  are  discussed.  Examples  are  presented  to 
examine  the  accuracy  of  these  proposed  approximations. 

In  Chapter  5,  we  describe  experiments  that  compare 
the  error  behaviour  of  various  approximations.  The  ranges 
of  the  parameter  values  and  the  configurations  of  the 
queueing  networks  are  explicitly  defined  with  the  help  of 
theoretical  results  and  actual  models.  We  shall  show  that 
one  of  the  causes  for  the  error  behaviour  in  some 


approximate  solution  techniques  is  due  to  the  particular 
representation  of  highly  skewed  service  time  distributions. 
This  chapter  serves  as  a  guideline  to  further  studies  into 
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the  behaviour  of  approximate  solution  techniques. 

Finally,  in  Chapter  6,  the  results  of  this 
are  summarized,  implications  of  this  work  on 
approximation  studies  are  discussed,  and  further 
directions  are  suggested. 


thesis 
f  ur  ther 
research 
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Chapter  2 

QUEUEING  NETWORK  MODELS  . 

In  this  thesis  we  are  concerned  with  the 
applicability  of  queueing  network  models  to  the  study  of  the 
performance  of  computer  systems.  The  basic  problems  we  must 
address  are  when  and  how  efficient  solution  techniques  can 
be  applied.  In  this  chapter  we  present  an  overview  of 
queueing  network  models  and  discuss  exact  and  approximate 
efficient  solution  techniques  that  can  be  used  to  solve 
them . 

The  first  section  introduces  notation  used  in 
queueing  theory,  and  illustrates  relationships  to  models  of 
computer  systems.  The  second  section  discusses  one  of  the 
most  important  assumptions  generally  needed  to  solve  a 
queueing  network  model  analytically,  namely,  the  exponential 
assumption.  The  third  section  discusses  both  recent 
advances  in  queueing  theory  and  efficient  exact  solution 
techniques  of  queueing  network  models.  Efficient 
approximate  solution  techniques  are  presented  in  the  fourth 
section.  This  section  is  divided  into  two  parts.  The  first 
part  discusses  approximations  based  on  aggregations  and  the 
second  part  deals  with  product  form  approximations. 
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2.1  Preliminaries 


A  queueing  network  is  composed  of  service  centres, 
each  consisting  of  a  queue  and  one  or  more  servers.  Each 
service  centre  is  connected  to  other  centres  by  transition 
paths.  Customers  require  service  at  service  centres.  The 
routing  of  the  customers  among  the  service  centres  is 
governed  by  transition  probabilities. 

A  service  centre  can  be  specified  by  an  arrival 
process,  a  service  time  distribution,  and  a  scheduling 
discipline  (also  known  as  service  discipline). 
Notationally ,  a  service  centre  is  represented  by  a/b/c/d, 
where  a  is  the  arrival  process,  b  is  the  service  process,  c 
is  the  number  of  parallel  servers  at  the  service  centre,  and 
d  is  the  maximum  number  of  customers  (requests)  that  can  be 
present  at  the  centre.  For  example,  an  M/M/1  centre  has 
Markovian  arrivals  (M) ,  a  Markovian  service  process  (M) ,  a 
single  server,  and  an  infinite  capacity  queue.  (When  the 
fourth  parameter  is  omitted,  it  is  assumed  to  be  infinity.) 
Markovian  processes  have  exponentially  distributed 
inter-event  times.  Similarly,  an  M/G/l/N  centre  has 
Markovian  arrivals  (M) ,  a  general  service  process  (G) ,  a 
single  server,  and  at  most  N  customers  present.  A  good 
discussion  of  various  single  queue  systems  can  be  found  in 
[Cooper  1973]. 

The  scheduling  discipline  determines  which 
customer (s)  in  the  queue  receive  service.  Some  queueing 
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disciplines  commonly  used  in  computer  systems  are:  first- 
come-f irst-served  (FCFS) ,  round-robin  (RR) ,  and  priority. 
For  example,  at  the  CPU,  a  scheduling  discipline  that 
attempts  to  give  fair  service  to  all  customers  is  RR.  Other 
disciplines,  which  are  used  in  queueing  network  models  for 
tr ac tabil ity ,  include 

last-come-f ir st-served-pr eemptive-r esume  (LCFSPR) ,  processor 
sharing  (PS)  (the  limiting  case  of  RR  when  the  quantum  size 
approaches  zero),  and  no-queueing  (NQ) .  In  NQ,  the  number 
of  servers  at  the  centre  is  greater  than  or  equal  to  the 
maximum  number  of  customers  that  can  be  at  the  service 
centre;  this  is  useful  for  modeling  service  at  a  group  of 
interactive  terminals. 

Customers  are  distinguishable  from  one  another  by 
the  amount  of  service  they  require  at  various  centres  and 
the  routing  pattern  they  follow  in  the  network.  Customers 
that  have  statistically  identical  behaviour  are  said  to 
belong  to  the  same  class.  A  queueing  system  that  has  more 
than  one  class  of  customers  is  termed  a  multiclass  system. 
For  example,  since  batch  and  interactive  jobs  almost  always 
have  different  resource  requirements,  they  are  typically 
classified  as  belonging  to  two  different  classes.  In  a 
network,  as  customers  move  from  one  centre  to  another,  they 
may  change  their  class  of  affiliation;  such  networks  are 
called  queueing  networks  with  class  changes.  The  amount  of 
service  time  a  customer  requires  at  a  service  centre  may 
depend  upon  its  customer  class  and  the  local  queue  length  at 
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the  centre.  Such  service  centres  are  called  class  dependent 
and  load  dependent  centres,  respectively. 

A  network  is  open  for  a  class  of  customers  if 
external  arrivals  and  departures  are  permitted  in  that 
class;  it  is  closed  for  a  class  of  customers  if  no  arrivals 
and  departures  are  allowed  in  that  class.  A  network  may  be 
open  for  some  classes  of  customers  and  closed  for  others 


such 

networks  are 

called  mixed 

.  Per 

haps  the 

most 

influential  queueing 

network  model 

is 

the 

central 

server 

model 

[Buzen  1971 ] , 

an  example  of 

wh 

ich  is 

shown  in 

Figure 

2. 1. 


I/O  Devices 


Figure  2.1:  Central  Server  Model 


2. 2  The  Exponential  Assumption 


Various  solution 
queueing  network  models 
combinations  of  the  two. 


techniques  can  be  used  to  solve 
analytic,  simulation,  and  hybrid 
For  the  reasons  mentioned  in 
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Chapter  1,  we  are  mainly  interested  in  efficient  analytic 
solution  techniques. 

A  queueing  network  model  can  be  characterized  by 
its  state.  The  set  of  all  feasible  states  is  known  as  the 
state  space.  Equilibrium  performance  measures  such  as 
server  utilization,  mean  queue  length,  queue  length 
distribution,  and  mean  response  time  can  easily  be 
calculated  from  the  equilibrium  state  probability 
distribution.  The  equilibrium  probability  distribution  can 
be  derived  by  solving  a  set  of  linear  equations  known  as  the 
global  balance  equations.  These  equations  relate  the  rates 
of  flow  into  and  out  of  each  state.  The  equilibrium 
assumption  sets  the  rate  of  flow  into  a  state  equal  to  the 
rate  of  flow  out  of  that  state.  For  any  network,  the  number 
of  global  balance  equations  that  must  be  solved 
simultaneously  is  the  same  as  the  size  of  the  state  space. 

Let  us  consider  the  M/M/1  system.  Because  of  the 
memoryless  property  of  the  exponential  distribution 
[Kleinrock  1975],  the  state  of  the  system  is  completely 
specified  by  the  number  of  customers  present. 

Exponential  servers  can  serve  as  building  blocks  in 
a  queueing  network  model.  In  a  single  class  network,  if  all 
the  servers  are  exponential,  the  state  of  the  network  is 
completely  described  by  the  number  of  customers  at  each 
server.  Let  us  consider  a  closed  network  with  N  jobs  and  M 
servers.  The  size  of  the  state  space  is  given  by 


The  state  space  therefore  grows  exponentially  with  M  and  N. 
For  example,  for  N=10  and  M=10,  the  number  of  states  in  the 
state  space  is  more  than  92,000  [Denning  and  Buzen  1978]. 

Two  basic  problems  need  to  be  addressed  at  this 

point : 

1.  The  exponential  assumption  is  too  strict;  that  is, 
observed  service  time  distributions  are  often  not 
exponentials.  More  importantly,  the  observed  distributions 
cannot  even  be  closely  approximated  by  exponentials. 

2.  The  number  of  equations  to  be  solved  grows  exponentially 
with  the  number  of  servers  and  customers. 

In  the  rest  of  this  section  we  discuss  solutions  to  the 
first  problem.  The  second  problem  is  addressed  in  the  next 
section. 

If  we  do  not  assume  an  exponential  service  time 
distribution  (e.g.,  consider  an  M/G/l  system),  the  state  of 
the  system  cannot  be  characterized  solely  by  the  number  of 
customers  at  each  centre  in  the  system.  The  state  also 
depends  upon  the  amount  of  service  already  received  by  the 
customer (s)  in  service;  it  cannot  be  solved  with  a  finite 
number  of  equations.  Thus,  the  global  balance  technique 
cannot  be  used  in  these  cases. 

A  method  generally  adopted  to  overcome  this  problem 
is  known  as  the  method  of  stages.  Any  c 
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distribution  with  a  rational  Laplace  transform  can  be 
represented  exactly  as  a  combination  of  exponential  servers 
[Cox  1955],  The  general  Cox-type  server  is  shown  in  Figure 

2.  2. 


Figure  2.2:  Cox-type  Server 


In  computer  system  modeling,  general  service  time 
distributions  have  usually  been  characterized  by  the  first 
two  moments  (however,  see  some  work  by  Lazowska  [1977]  and 
Bux  and  Herzog  [1977]).  The  coefficient  of  variation  (CV) 
of  a  distribution  is  defined  as  the  variance  divided  by  the 
square  of  the  mean.  Given  any  mean  (M)  and  positive  CV,  a 
distribution  with  the  specified  first  two  moments  can  be 
found  with  one  of  the  following  three  forms  [Sauer  and 
Chandy  1975,  Sevcik  et  al.  1977]  (Figure  2.3): 

HE  2 :  Hyperexponential  (2-stage)  (CV>1) 
f (t) =a*r*exp(-r*t)  +  (1-a) *s*exp (-s*t) 

where  a=  (1/2) * [ 1 -sqr t (  (CV-1 ) / (CV+1 ) ) ] , 
s=2* ( 1-a ) /M ,  and 
r  =  2  *a/M 
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Exponential 


Generalized  Erlang  (k-stages) 


Figure  2.3:  Distributional  Forms 
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EXP:  Exponential  (CV=1) 
f ( t) =r*exp (-r *t) 
where  r=l/M 

GEK :  Generalized  Erlang  (K-stage)  (CV<1) 

K-i 

f (t)=(r*t) (l-a+(a*r*t/(K-l) ) ) *r*exp (-r*t) / (K-2) ! 
where  K=.fl/CVl, 

T=l- ( K  —  1 ) *CV 
a=T+sqrt(K*T)/(CV+l) , 
r  =  (K-l+a ) /M ,  and 

The  size  of  the  state  space  grows  exponentially 
with  the  number  of  stages  one  has  to  use  to  represent  a 
service  time  distribution  at  a  service  centre.  Thus,  the 
number  of  equations  to  be  solved  grows  exponentially  with 
the  number  of  stages  used  to  represent  service  time 
distributions.  In  certain  cases,  however,  we  can  get  around 
the  problem  of  solving  a  large  number  of  equations  and 
calculate  exact  solutions  efficiently.  The  assumptions  one 
has  to  make  to  obtain  exact  efficient  solutions  are 
discussed  in  the  next  section. 

2.3  Efficient  Exact  Solution  Techniques 


Jackson  [1957,  1963]  studied  mixed,  single  class 
networks  composed  of  exponential  service  centres  at  which 
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service  rates  may  depend  upon  the  local  queue  length. 
Because  Jackson  considered  only  exponential  servers,  a  state 
of  the  system  could  be  represented  as  (nl,n2, . . . ,nM) ,  where 
ni  is  the  number  of  customers  present  at  the  i-th  service 
centre  and  M  is  the  total  number  of  centres.  He  showed  that 
the  equilibrium  solution  to  such  a  network  has  a  product 
form.  That  is,  the  equilibrium  solution  for  the  probability 
that  the  system  is  in  a  given  state  is  the  product  of  the 
independent  probabilities  of  the  individual  servers  being  in 
their  respective  states: 


P  [ n 1 , n  2 , . . . , n M ]  =  (1/G) *P [nl ] *P [n2] *. . . *P  [nM] 


where  G  is  a  normalization  constant 
probabilities  for  all  the  states 
Newell  [1967]  developed  a  separation 
technique  for  closed  networks  wi 
times.  P(ni)  is  calculated  based  on 
rate  to  the  ith  server.  If  (rl,r2, 
arrival  rates  then  they  must  satisfy 
equations  : 

M 

ri  Pji*rj  for 

j  =  l 

where  Pij  is  the  branching  probabili 
The  above  set  of  equations,  als 
constraints,  relate  the  throughputs 
The  normalization  constant, 
relative  probabilities  for  all  the 
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and  Kobayashi  1976].  Relative  probabilities  are  calculated 
based  on  the  relative  arrival  rates  to  individual  servers. 
In  open  systems  absolute  arrival  rates  are  known  precisely. 
This  reduces  computational  problems  with  respect  to  G;  in 


such  systems  G=l. 

In  general,  open  networks  require  no 

computation  for 

G.  G  is  meaningful  only  for  a  closed 

network  where  the 

total  number  of  customers  is  fixed,  say  N, 

and  G  is  represented  as  a  function  of  the  fixed  number  of 
customers  N,  G(N). 


A  naive 

approach  to  calculating  G(N)  is  to 

enumerate  all 

states  and  calculate  their  relative 

probabilities . 

Absolute  probabilities  can  then  be 

determined  from 

the  relative  probabilities  by  normalizing 

their  sum  to  one. 

It  is  important  to  note  that  the  rapid 

growth  of  the  state  space  makes  this  infeasible  for  all  but 
the  smallest  networks. 

Buzen  [1971,  1973]  presented  an  efficient  algorithm 
to  calculate  G(N)  without  enumerating  all  the  states  of  the 
network.  His  convolution  approach  to  evaluating  G(N)  is 
very  efficient  and  requires  little  storage.  For  example, 
total  storage  requirement  is  2*N  locations,  and  the 
algorithm  has  an  operations  count  of  about  2*M*N  if  all  the 
servers  are  load  independent,  and  2*M*N*(N+1)  otherwise  -  a 
large  saving  compared  to  enumerating  all  the  states. 


The  exponential  assumption  poses  a  major  drawback 
to  the  application  of  queueing  network  models  to  computer 
systems.  Service  time  distributions  at  various  servers  in 
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computer  systems  are  rarely  exponential  [Baskett  1971, 
Lester  1973].  In  particular,  the  service  time  distributions 
are  often  found  to  have  a  high  variance  (CV>>1)  at  the  CPU 
and  a  low  variance  (CV<1)  at  the  I/O  devices.  Assuming 
exponential  service  time  distributions  in  such  cases  can 
introduce  large  errors  in  the  calculations  of  performance 
measures  [Chandy  et  al.  1975b]. 

The  assumption  of  exponential  service  time 
distributions,  however,  was  basic  to  Jackson's  product  form 
solution,  which  in  turn  was  essential  to  making  exact 
solutions  of  such  systems  tractable.  Baskett  [1971]  showed 
that  the  equilibrium  distribution  of  system  state  of  a 
network  with  servers  having  service  time  distributions  with 
rational  Laplace  transforms  and  PS  scheduling  disciplines  is 
the  same  as  that  for  a  network  where  the  servers  have 
exponential  service  time  distributions.  (Note  that  for  a 
single  class  of  customers  with  exponential  service  time 
distributions,  as  long  as  a  server  is  not  idle  when  there  is 
a  customer  present  at  the  service  centre,  the  scheduling 
discipline  has  no  effect  on  the  equilibrium  distribution  of 
system  state.) 

Chandy  [1972]  proved  that  networks  with  service 
time  distributions  with  rational  Laplace  transforms  and  PS 
or  LCFSPR  service  disciplines  also  have  product  form 
solutions.  These  solutions  depend  only  upon  the  mean  of  the 
service  time  distributions  and  not  upon  the  higher  moments. 
Chandy's  results  were  derived  using  the  property  known  as 
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local  balance.  Informally,  a  network  is  in  local  balance  if 
the  rate  of  flow  out  of  a  state  due  to  a  departure  from  a 
specific  service  centre  is  eaual  to  the  rate  of  flow  into 
that  state  due  to  an  arrival  to  the  same  service  centre  for 
each  class. 

Baskett  et  al .  [1975]  unified  these  results  and 
showed  that  the  concept  of  local  balance  could  be  extended 
to  multiple  classes  and  other  types  of  scheduling 
disciplines.  Thus,  a  large  class  of  networks  that  have 
product  form  solutions  was  established.  A  network  has 
product  form  if  each  service  centre  is  one  of  the  following 
types : 


1.  FCFS  where  all  customers  have  the  same  exponential 
service  time  distribution;  servers  may  be  load  dependent, 

2.  PS  and  a  service  time  distribution  with  rational  Laplace 
transform , 


3.  LCFSPP  and  a  service  time  distribution  with  rational 
Laplace  transform,  or 

4.  NQ  and  a  service  time  distribution  with  rational  Laplace 
transform. 

The  product  form  solution  for  multiple  classes  is  given  by 
P  [SI, S2, . . . , SM]  =  (1/G (N) )*P[S1]*P[S2]*...*P [SM] 


where  Si  =  (nil,ni2, . . . ,niK) , 

nij  is  number  of  class  j  customers  at  centre  i,  and 
K  is  the  number  of  classes  in  the  network. 


The  results  derived  by  Baskett  et  al.  [1975]  allow 
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class  changes,  but  assume  that  the  routing  of  the  jobs 
follows  a  first-order  Markov  chain  (i.e.,  transitions  are 
memoryless).  Reiser  and  Kobayashi  [1975]  allow 
transitions  among  service  centres  to  be  governed  by  an 
arbitrary  (finite)  order  Markov  chain.  By  using  artificial 
class  changes  in  conjunction  with  a  single  order  Markov 
chain  to  retain  more  information  about  the  past  routing  of 
the  jobs,  a  higher  order  Markov  chain  can  be  viewed  as  a 
single  order  chain. 

Various  other  sets  of  sufficient  conditions  for  a 
system  to  have  a  product  form  solution  have  also  been 
derived.  We  mention  some  of  these  here  briefly. 


1.  Station  balance  and  local  balance:  Chandy  et  al.  [1977] 
define  a  queue  as  a  set  of  stations,  which  can  be  occupied 
by  customers  present  at  that  centre.  A  queueing  discipline 
satisfies  station  balance  if  the  rate  at  which  a  customer 
receives  service  at  a  station  of  the  queue  is  proportional 
to  the  probability  that  the  customer  arrives  to  that 

station.  Station  balance  is  a  local  property  of  a  queue  at 
the  service  centre  whereas  local  balance  is  an  overall 
property  of  a  service  centre.  It  is  stronger  than  local 
balance  in  the  sense  that  it  imposes  restrictions  on  each 
position  of  the  queue.  Chandy  et  al .  [1977]  show  that 
station  balance  at  each  server  in  a  network  is  sufficient  to 
guarantee  that  the  network  has  a  product  form  solution. 


2.  M  =>  M  property: 


If  a  Markovian  arrival  process  to  a 
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server  implies  a  Markovian  departure  process, 
said  to  have  M  =>  M  property.  Muntz  [1973] 
networks  with  servers  satisfying  the  M  =>  M 
product  form  solutions  and  that  the  networ 
satisfy  the  M  =>  M  property- 


the  server  is 
showed  that 
property  have 
ks  themselves 


3.  Lost  and  triggered  arrivals:  Lam  [1977]  has  shown  that 
under  certain  conditions  a  queueing  network  with  state 
dependent  external  arrivals  can  also  have  a  product  form 
solution.  Depending  on  the  state  of  the  system  (i.e., 
either  total  number  of  customers  in  the  system  or  numbers  in 
a  particular  closed  chain)  an  external  arrival  enters  the 
system  or  is  lost.  Also,  for  certain  states  a  departure 
from  the  system  causes  immediate  arrival  (triggered)  into 
the  system  of  another  customer  of  the  same  class.  Lam  gives 
sufficient  conditions  for  such  a  system  to  have  local 
balance  (and  therefore  product  form). 

Various  sets  of  sufficient  conditions  for  a  system 
to  have  a  product  form  have  some  relationships  among 
themselves.  Figure  2.4  presents  known  implications  [Chandy 
et  al.  1977]  among  some  of  these  conditions. 

The  convolution  algorithm  mentioned  earlier  has 
been  extended  to  include  the  networks  described  in  [Baskett 
et  al.  1975]  by  Reiser  and  Kobayashi  [1975],  Wong  [1975], 
and  Chandy  et  al.  [1975a].  This  algorithm  has  been 
generalized  to  include  a  more  general  class  of  networks  with 
lost  and  triggered  arrivals  by  Zahorjan  [1979].  Good 


discussions  of  the  algorithm  to  calculate  the  equilibrium 
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Figure  2.4 
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probabilities  can  be  found  in  [Shum  1976,  Reiser  and 
Kobayashi  1976,  and  Reiser  and  Sauer  1977].  A  number  of 
solution  packages  using  this  technique  have  been  developed 
(e.g.,  ASQ  [Keller  1973],  BEST/1  [Buzen  1978],  QNET-4 
[Reiser  and  Sauer  1978]). 

Unfortunately,  there  are  some  very  basic  types  of 
service  centres  in  computer  system  models  that  do  not 
satisfy  the  conditions  for  product  form  solutions. 
Moreover,  these  types  of  service  centres  are  observed  in 
many  actual  systems.  For  example,  a  FCFS  scheduling 
discipline  with  a  non-exponential  service  time  distribution, 
which  we  observe  in  many  actual  systems,  does  not  fall  into 

this  category,  even  for  a  single  class.  In  the  next  section 
we  present  several  approximation  techniques  that  give 

efficient  solutions  to  networks  that  do  not  satisfy  the 
conditions  for  a  product  form  solution. 


2.4  Efficient  Approximate  Solution  Techniques 


If  a  queueing  network  does  not  satisfy  conditions 
for  a  product  form  solution,  we  cannot  use  the  efficient 
solutions  discussed  in  the  previous  section  to  get  an  exact 
solution.  In  these  situations,  to  get  an  efficient 
solution,  we  are  left  with  two  alternatives: 

1.  Assume  the  conditions  required  for  product  form  and  then 
solve  the  model  using  an  efficient  exact  solution  technique. 
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Although  the  solution  technique  is  exact,  we  are  solving  a 
model  that  approximates  the  original  model. 

2.  Solve  the  model  using  an  approximate  solution  technique. 

There  are  two  basic  disadvantages  to  the  first 
approach:  (1)  it  can  produce  very  large  errors  in  the 
performance  measures  [Chandy  et  al.  1975b],  and  (2)  the 
necessary  assumptions  limit  the  predictive  ability  of  the 
model.  For  example,  if  the  CPU  is  FCFS  with  a 
non-exponential  service  time  distribution,  we  must  assume  an 
exponential  distribution  (or  a  PS  scheduling  discipline)  for 
a  product  form  solution.  This  eliminates  the  possibility  of 
studying  the  behaviour  of  system  performance  with  respect  to 
large  variations  in  CPU  service  time  distributions. 

Approximate  solution  techniques  also  have  the  first 
disadvantage  mentioned  above,  but  are  free  from  the  second 
one.  Although  we  still  have  the  possibility  of  large 
errors,  we  can  attempt  to  improve  an  approximation  to 
minimize  possible  errors,  unlike  the  first  case.  Chapter  3 
addresses  the  issue  of  improving  approximate  solution 
techn iques . 

There  are  two  basic  approaches  to  calculating  an 
efficient  approximate  solution  of  a  queueing  network  that 
does  not  have  product  form  solution: 

1.  Aggregate  parts  of  the  network  and  reflect  the  combined 
behaviour  by  a  single  server,  thus  reducing  the  complexity 
of  the  network.  This  is  an  approximation  because  the 
single  server  may  not  accurately  reflect  the  part  of  the 
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network  it  is  replacing-  Such  aggregations  could  be  applied 
repeatedly. 

2-  Modify  the  solutions  of  the  individual  servers  to  reflect 
the  face  that  the  conditions  for  product  form  solution  are 
not  satisfied,  but  still  use  product  form-  This  is 
approximate  because  the  parameters  to  solve  individual 
service  centres  may  depend  upon  the  behaviour  of  the  rest  of 
the  network-  We  will  call  these  product  form 
approximations- 

The  first  approach  solves  multiple  centres  together 
and  the  second  approach  solves  individual  centres,  before 
combining  these  individual  results  to  get  the  overall 
solution-  In  the  remainder  of  this  section,  both  these 
approaches  are  discussed  with  the  help  of  some  commonly  used 
approximations,  and  some  other  approximations,  which  do  not 
belong  to  these  two  groups,  are  presented. 

2-4.1  Aggregation/Decomposition 

Aggregation  and  decomposition  are  simply  the 
complementary  intuitive  processes  of  combining  and 
separating  parts  of  the  system  to  facilitate  analysis. 
Aggregation  involves  solving  portions  of  the  model  in 
isolation  and  gathering  results  together  to  obtain  a 
solution  of  the  entire  model-  After  analyzing  the 
subnetworks,  they  are  each  replaced  by  a  single  composite 
server,  forming  the  reduced  model.  A  composite  server  is  a 
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single  load  dependent  server  that  reflects  the 
characteristics  of  the  subnetwork  it  replaces.  In  this  way, 
the  complexity  of  solving  the  network  is  reduced 
considerably.  The  process  of  dividing  a  system  into  many 
parts,  where  each  part  can  be  solved  separately,  is  known  as 
decomposition.  In  solving  each  part  of  the  system,  the 
behaviour  of  that  part  of  the  system  is  represented  by  a 
single  aggregated  (composite)  server.  As  mentioned  earlier, 
this  can  introduce  error  because  it  may  not  be  possible  to 
reflect  the  exact  behaviour  of  the  subsystem  by  a  single 
server  with  exponentially  distributed  inter -depar ture  times. 
In  addition  to  being  an  approximation,  aggregation  can  also 
limit  the  predictive  ability  of  the  model  since  performance 
measures  for  individual  service  centres  in  the  subnetwork 
may  not  be  available  from  the  overall  solution  of  the 
network.  Recently,  Sauer  [1979]  has  presented  a  method  for 
evaluating  the  performance  measures  for  individual  centres 
in  aggregated  subsystems.  The  method  gives  exact  solutions 
for  the  networks  that  satisfy  the  conditions  for  a  product 
form.  Applicability  of  these  results  in  other  networks 
needs  further  investigation. 

The  simplification  provided  by  decomposition  may  be 
more  than  simply  useful:  it  may  be  necessary.  Even  for 
systems  satisfying  conditions  for  a  product  form  solution, 
as  the  network  becomes  larger  many  solution  techniques 
experience  numerical  instabilities.  For  example,  when  there 
are  many  jobs,  classes,  and  load  dependent  devices,  the 
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numbers  computed  to  obtain 
overflow  even  with  double 
Buzen  1978]. 


the  normalization  constant  may 
precision  arithmetic  [Denning  and 


DECOMPOSITION 


Intuitively,  a  system  can  be  decomposed  into 
smaller  subsystems  if  (i)  interactions  within  subsystems  can 
be  studied  as  if  interactions  among  subsystems  did  not  exist 
and,  (ii)  interactions  among  subsystems  can  be  analyzed 
without  referring  to  the  interactions  within  subsystems 
[Courtois  1975,  1977].  A  system  for  which  these  conditions 
hold  is  said  to  be  nearly-decomposable . 

To  illustrate  near-decomposabil ity  let  us  consider 
a  closed  queueing  system  composed  of  two  subsystems  SI  and 
S2  (Figure  2.5).  Jobs  completing  service  at  a  subsystem  are 
fed  back  to  it  with  probability  P  and  branch  to  the  other 
subsystem  with  probability  1-P.  As  P  approaches  1,  the  rate 
of  flow  of  customers  between  the  two  subsystems  approaches 
zero.  A  system  with  subsystems  that  have  very  small  rates 
of  flow  of  jobs  among  the  subsystems  is  called  weakly 
coupled . 


The 

solution  for 
by  assuming 
inter-subsys 


effort  in  calculating  the  overall  equilibrium 
a  weakly  coupled  system  is  reduced  considerably 
that  the  subsystems  are  in  equilibrium  between 
tern  interactions.  The  service  time  distribution 


(inter-departure  times)  at  each  subsystem  need  only  to  be 
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calculated  for  various  numbers  of  jobs  present  in  the 
subsystem. 

The  service  time  distribution  of  a  subsystem  will 
be  the  inter-departure  distribution  of  that  subsystem.  As 
the  rate  of  flow  of  customers  between  subsystems  approaches 
zero,  the  inter-departure  distribution  approaches 
exponential.  Thus,  in  estimating  the  service  time 
distribution  of  the  subsystem,  we  need  only  estimate  the 
mean,  not  the  higher  moments  of,  inter-departure  times  for 
the  subsystem. 

If  a  subsystem  interacts  infrequently  with  its 
environment,  the  transient  behaviour  of  the  subsystem  will 
have  little  effect  on  the  long  term  (equilibrium)  behaviour 
of  the  total  system.  Thus,  very  little  error  will  be 
introduced  by  assuming  that  the  subsystem  is  in  equilibrium 
for  the  entire  interval  between  two  interactions  with  the 
environment. 

Simon  and  Ando  [1961]  and  Ando  and  Fisher  [1963] 
proved  that  for  a  given  standard  of  accuracy  there  always 
exists  a  sufficient  degree  of  near-decomposabil ity  such  that 
an  aggregation  procedure  will  meet  the  standard.  It  should 
be  observed  that  for  a  high  standard  of  accuracy  or  in  a 
system  that  is  strongly  (as  opposed  to  weakly)  coupled  the 
amount  of  reduction  in  the  size  of  the  overall  model  may  be 
minimal . 

Courtois  [1975,  1977]  addressed  the  problem  of  how 
much  error  is  incurred  by  decomposing  a  system  in  a  given 
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way.  As  one  would  expect,  the  error  in  using  decomposition 
increases  with  the  rate  of  interactions  among  the 
subsystems.  Since  Courtois  dealt  with  solutions  of  linear 
equations  (global  balance  equations) ,  well-developed  error 
analysis  techniques  in  numerical  analysis  can  be  used  to 
estimate  the  error  introduced  by  decomposition.  Courtois 
[1977]  presents  a  set  of  necessary  and  sufficient  conditions 
for  near-decomposability  and  a  technique  to  estimate  the 
error  introduced  due  to  decomposition.  This  is  one  of  the 
most  important  advantages  of  this  approximation  over  others. 
Zarling  [1976],  however,  shows  that  the  constant  of 
proportionality  for  the  error  bounds  may  be  large  enough  to 
prohibit  use  of  the  technique  in  many  cases.  The  error  is 
also  found  to  depend  in  a  critical  way  on  the  particular 
choice  of  how  to  decompose  the  state  transition  matrix 
describing  the  system.  The  usefulness  of  the  technique  is 
limited  by  the  size  of  the  state  transition  matrix,  which  in 
general  is  large.  Also,  it  is  unclear  whether  Courtois' 
numbers  are  either  good  estimates  or  firm  bounds  [Lazowska 
and  Garner  1979 ] . 

As  an  example  of  a  computer  system  where  near- 
decomposability  can  be  used,  consider  an  interactive 
computer  system  (Figure  2.6).  The  terminals  can  be  viewed 
as  one  service  centre  with  NQ.  The  interactions  between  the 
terminals  and  the  system  are  rare  compared  to  interactions 
among  servers  inside  the  system.  Thus  we  can  decompose  the 
interactive  system  into  two  subsystems  -  the  terminals  and 
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Figure  2.5:  Near-Decomposability  of  a  System 


Terminals 


Figure  2.6:  Interactive  System 
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the  central  subsystem. 

We  conclude  by  noting  that  near -decomposabil i ty 
only  provides  a  good  approximation  in  certain  circumstances. 
For  example,  the  central  subsystem  without  the  terminals  may 
itself  be  too  large  to  solve  using  global  balance 
techniques.  If  the  interactions  between  the  CPU  and  I/O 
devices  are  very  frequent,  assuming  near -decomposabil ity  of 
the  CPU  and  I/O  devices  can  produce  large  errors. 

NORTON  THEOREM  REDUCTION 

Chandy  et  al.  [1975a]  proposed  an  aggregation 
technique  in  which  a  portion  of  a  queueing  network  can  be 
"short-circuited"  and  considered  in  isolation  as  a  closed 
network.  From  this  we  can  calculate  the  load-dependent 
throughput  of  the  closed  system,  which  defines  the  service 
time  (inter-departure  time)  of  the  composite  server.  For 
example,  the  queueing  network  model  represented  in  Figure 
2.7a  reduces  to  the  model  given  in  Figure  2.7b  by  replacing 
the  subsystem  in  the  box  by  a  composite  server. 

Replacing  part  of  the  network  by  a  server  with  only 
the  mean  output  rate  specified,  however,  is  the  same  as 
assuming  an  exponential  inter-departure  time  distribution  at 
the  composite  server.  This  approach  can  introduce  a  large 
error  in  the  solution  of  the  system. 

Sauer  and  Chandy  [1975]  suggested  a  heuristic  for 
estimating  the  service  time  distribution  of  the  composite 
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server  more  accurately  than  by  just  assuming  an  exponential 
distribution.  They  proposed  that  for  parallel  servers  the 
CV  of  the  composite  server  be  the  weighted  average  of  the 
CV's  (where  the  weights  are  the  branching  probabilities)  of 
the  component  servers.  For  example,  consider  the  subsystem 
shown  in  Figure  2.8.  Let  Pi  be  the  branching  probability  to 
the  i-th  server,  which  has  a  service  time  distribution 
having  mean  Mi  and  coefficient  of  variation  CVi.  The  CV  of 

the  composite  server  (CVc)  is  given  by: 

M 

CVc  =y^Pi*CVi  (2.1) 

L*l 

Then,  depending  upon  CVc,  a  distributional  form  is  chosen 
for  the  composite  server. 

This  approximation  is  in  some  ways  intuitively 
pleasing,  and  in  many  situations  gives  satisfactory  results. 
It  is  not  difficult,  however,  to  identify  systems  for  which 
the  approximation  leads  to  high  errors  (more  than  20%)  in 
server  utilizations  [Tripathi  and  Levy  1976]. 

Strict  modeling  of  the  aggregate  service  time 
produces  the  following  formula 

H  2  M  z  M  m 

CVc  =  [£Pi(CVi*Mi  )  +£pi*Mi  -  ]T(Pi*Mi  )2  ]  /  (£Pi*Mi  )z 

L~l  c-i  »'=/ 

This  formula  yields  a  coefficient  of  variation  that  is 
exactly  the  coefficient  of  variation  of  the  overall 
distribution  of  service  times  in  the  subsystem  [Tripathi 
1976].  Experiments  have  shown  that  using  this  exact  formula 
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Figure  2.7a 


Figure  2.8:  I/O  System  Aggregation 
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for  CVc  can  still  lead  to  unacceptable  error  and  thus  it 
does  not  give  a  uniform  improvement  over  the  one  suggested 
by  Sauer  and  Chandy  [Tripathi  and  Levy  1976]. 

Sevcik  et  al.  [1977]  suggested  that  it  may  not  be 
enough  to  improve  the  accuracy  of  solving  the  composite 
server  locally  without  any  consideration  of  interaction  with 
the  rest  of  the  system.  The  composite  server  should  reflect 
the  same  behaviour  externally  as  did  the  subsystem.  In 
other  words,  given  an  arrival  process  to  the  subsystem,  the 
composite  server  should  approximate  as  closely  as  possible 
the  departure  process  of  the  original  subsystem.  This  is 
referred  as  the  input/output  characterization  issue. 

In  Chapter  3,  we  show  that  this  is  a  common  problem 
with  most  of  the  approximate  solution  techniques.  To 
characterize  the  input/output  event  stream,  we  must  address 
three  subproblems  (see  Figure  2.9): 

1.  How  does  splitting  affect  an  event  stream? 

2.  What  is  the  effect  on  an  event  stream  of  going  through  a 
server? 

3.  What  is  the  result  of  merging  event  streams  together? 

The  last  two  questions  are  very  hard  to  answer,  although 
Sevcik  et  al.  [1977]  and  Kuhn  [1976]  provide  approximate 


solutions . 
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Figure  2.9:  Input/Output  Characterization 
Approximate  solution  techniques  have  also  been 
proposed  for  queueing  network  models  with  priority 
scheduling  disciplines.  Sauer  and  Chandy  [1975]  reduce  the 
complexity  of  the  model  by  aggregating  classes  and 

aggregating  I/O  subsystems  (using  Norton's  theorem  for 

queueing  networks).  After  aggregation  of  classes  and 

centres r  the  reduced  model  is  sufficiently  small  to  be 
solved  quickly  even  though  it  is  not  in  local  balance. 


2.4.2  Product  Form  Approximations 


For  networks  that  have  product  form  solutions,  the 
equilibrium  probability  that  the  network  is  in  a  given  state 
is  the  product  of  individual  equilibrium  probabilities  for 
each  service  centre.  In  these  situations  the  second  and 
higher  moments  of  the  service  time  distribution  do  not  play 
any  part  in  the  solution.  The  individual  solutions  for  each 
service  centre  are  solutions  of  M/M/c/N  systems. 
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If  the  conditions  for  product  form  are  not 
satisfied,  the  solutions  of  the  individual  service  centres 
can  be  modified  to  obtain  an  efficient  approximate  solution 
of  the  network: 

P[nl,n2, ... ,nM]  =  ( 1/G (N) )*Fl(nl)*F2(n2)*...*FM(nM) 

where  Fi(ni)  is  the  modified  solution  (in  place  of  an 
M/M/c/N  solution)  for  the  .i-th  service  centre. 

There  have  been  three  basic  approaches  used  to 
obtain  Fi (ni ) : 

1.  Diffusion  approximations, 

2.  Analytic  solutions,  and 

3.  Numerical  approximations. 

In  the  remainder  of  this  section  these  three  approaches  are 
presented . 

DIFFUSION 

In  diffusion  approximations  every  server  is  viewed 
as  a  G/G/l  system.  G/G/l  systems  have  very  complex  analytic 
solutions,  which  are  expensive  to  evaluate.  Diffusion 
approximations  to  these  G/G/l  systems  are  used  for  Fi(ni). 
Once  solutions  to  the  G/G/l  systems  are  obtained,  a  product 
form  solution  algorithm  is  used  to  obtain  the  solution  for 
the  entire  network. 

In  a  diffusion  approximation,  queue  length  is 
treated  as  a  continuous  rather  than  discrete  value.  The 
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change  in  the  queue  length  in  a  small  interval  is 
approximated  by  a  Gaussian  process.  The  inter -ar r ival 
process  and  service  process  are  approximated  by  Gaussian 
processes  based  upon  their  first  two  moments,  the  mean  and 
the  variance.  In  general,  this  should  increase  the  accuracy 
over  assuming  exponential  inter-event  time  distributions  and 
taking  only  the  mean  into  consideration.  Newell  [1971] 
gives  an  excellent  introduction  to  and  motivation  for 
diffusion  approximations. 

There  are  two  basic  problems  in  using  diffusion 
approximations  for  a  G/G/l  system  -  boundary  conditions  and 
mapping  continuous  queue  lengths  onto  discrete  queue 
lengths.  We  must  impose  conditions  on  the  boundary  -  the 
empty  state  (no  customer  at  the  centre).  Different 
methods  have  been  proposed  for  this  problem.  Kobayashi 
suggested  using  a  reflecting  boundary,  i.e«,  whenever  the 
queue  hits  the  empty  state  (the  boundary),  an  instantaneous 
return  from  it  is  assumed  [Kobayashi  1974].  This  means  that 
the  server  is  assumed  to  be  heavily  utilized.  Gelenbe 
[1975,  1976]  has  shown  that  a  general  holding  time  at  the 
boundary  can  be  assumed,  and  that  only  the  first  moment  of 
the  holding  time  appears  in  the  solution.  The  probabilities 
obtained  by  a  diffusion  approximation  are  for  continuous 
queue  lengths.  Three  methods  for  estimating  the  probability 
for  queue  length  k,  namely,  sum  for  the  interval  (k-l,k), 
(k,k+l),  or  (k-1/2, k+1/2) ,  have  been  successfully  used  [Yu 
1977] . 
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The  calculations  involved  in  diffusion 
approximations  are  very  simple.  For  example,  let  Ma  and  CVa 
represent  the  mean  and  the  CV  of  the  inter-arrival  time,  and 
Ms  and  CVs  be  the  mean  and  the  CV  of  the  service  time. 
Also,  let  u=Ms/Ma  and  E  [q ]  be  the  mean  queue  length.  Then 


E[q]  =  u/(l-ur) 

where  u**  =  exp  (-2*  (1-u)  /  (CVs+u*CVa)  )  [Kobayashi  1974] 


or , 


E[q]  =  u* (l+(u*CVa+CVs)/2* (1-u) )  [Gelenbe  1975] 


The  application  of  diffusion  approximations  to 
general  queueing  networks  was  first  done  by  Kobayashi 
[1974].  The  state  of  the  network  is  defined  by  the  queue 
length  at  each  server.  The  interactions  among  different 
queueing  processes  are  explicitly  considered  in  the 
diffusion  equations  in  terms  of  the  second  moments  of  the 
arrival  processes  at  the  service  centres.  In  the 

resulting  solution,  the  queue  length  distribution  is 
expressed  as  a  product  of  marginal  queue  length 
distributions  at  each  server.  Thus  the  problem  is 
simplified  and  each  queue  can  be  solved  individually  by 
characterizing  the  input/output  process  accurately. 

Because  each  queue  is  analyzed  individually,  the 
accuracy  of  estimating  the  CV's  of  the  arrival  and  departure 
processes  at  each  server  plays  an  important  role  in  the 
accuracy  of  the  approximation.  Two  methods  for  estimating 


the  mean  arrival  rates  are  suggested  [Reiser  and  Kobayashi 
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1975]  :  assume  that  a  product  form  exists  and  solve  the 
network  to  get  the  arrival  rates,  or  assume  that  the  most 
heavily  utilized  server  is  always  busy  (its  utilization  is 
1).  (The  second  method  gives  upper  bounds  for  the  arrival 
rates.  The  first  method,  however,  does  not  provide  any  such 
bounds.  The  results  from  the  first  method  may  fall  above  or 
below  the  actual  arrival  rates  depending  upon  such 
parameters  as  the  CV's  of  the  service  time  distributions.) 
These  approximations  can  introduce  large  errors.  These 
errors  are  discussed  in  Chapters  3  and  4. 

Different  methods  of  estimating  these  CV's  have 
been  suggested.  Let  CVai,  CVsi,  and  CVdi  be  the  CV's  of 
the  inter-arrival  time,  the  service  time,  and  the 
inter-departure  time  at  the  ith  server,  respectively.  Also, 
let  Mai  be  the  mean  inter-arrival  time  at  the  ith  server  and 
Pij  be  the  transition  probability  from  server  i  to  j.  Then: 
1.  Kobayashi  [1974]  approximated  CVdi  as  CVsi,  i.e.,  a  fully 
utilized  server  is  assumed.  Then  CVai  is  represented  as 
follows : 


CVai  =  Mai*^[  (CVd  j-1)  *Pj  i  +  1  ]  *  (1/Ma  j  )  *Pji  (2,2) 

jxO 

where  M  is  the  total  number  of  servers  in  the  network 
(server  0  is  the  external  arrival  source  in  open  networks). 
2.  Gelenbe  and  Pujolle  [1976]  tried  to  remove  the  assumption 
that  servers  are  fully  utilized  and  gave  an  approximation  to 
CVdi  as 


CVdi 


,  .2 

1+ui  * 


(CVsi-1) + (1-ui) * (CVai-1) 


(2.3) 
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where  ui  is  the  utilization  of  the  i-th  server.  (As 
mentioned  before,  in  closed  networks  ui  can  be  determined 
either  by  assuming  that  the  network  has  a  product  form  or  by 
assuming  that  the  most  heavily  utilized  server  is  always 
busy.)  To  get  CVai  now,  we  have  to  solve  a  set  of  linear 
equations.  This  approximation  gives  better  results  than 
those  obtained  by  Kobayashi's  method  in  1  above. 

3.  Yu  has  suggested  a  "simpler"  approximation  to  CVdi  [Yu 
1977] : 

CVdi  =  1+ui2*  (CVsi-1) 

That  is,  an  M/G/l  queueing  system  at  each  server  is  assumed. 
Although  this  approximation  is  inexpensive  to  compute,  it 
ignores  the  arrival  processes  at  the  service  centres. 
Extensive  comparisons  have  not  yet  been  made  to  determine  if 
using  the  simple  form  still  maintains  reasonable  accuracy. 
We  discuss  this  problem  in  Chapter  3. 

Gelenbe  and  Pujolle  [1977]  have  extended  the 
applicability  of  diffusion  approximations  to  networks  with 
multiple  classes.  At  the  moment  this  analysis  is  limited  to 
open  networks  only.  In  Chapter  6,  one  way  of  extending  this 
analysis  to  closed  networks  is  mentioned. 

ANALYTIC 

Although  diffusion  approximations  are  broadly 
applicable,  the  large  number  of  approximations  required  to 
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get  a  G/G/l  solution  may  produce  a  large  error  in  the 
overall  solution.  Shum  and  Buzen  [1977]  assumed  that  the 
complement  of  each  server  is  a  composite  server  with 
exponentially  distributed  service  times,  and  that  therefore 
M/G/l  solutions  can  be  used  in  open  systems.  This  allows 
them  to  use  exact  techniques  to  analyze  each  server 
separately.  In  closed  systems,  M/G/l/N  solutions  are 
proposed  for  the  servers. 

Shum  [1976]  derives  expressions  for  M/G/l/N 
solutions  in  special  cases.  A  general  service  time 
distribution  is  classified  as  one  of  three  forms 
(exponential,  hyperexponential,  and  hypoexponential )  based 
upon  the  value  of  its  first  two  moments.  (The  method  has 
already  been  discussed  in  section  2.2.)  Closed  form 
expressions  for  these  three  cases  are  given  by  Shum. 

The  absolute  arrival  rate  at  each  server  is  needed 
before  M/G/l/N  solutions  can  be  obtained.  (Recall  that  this 
problem  in  not  difficult  in  open  systems.)  Shum[1976] 
suggests  an  iterative  technique  to  estimate  the  arrival 
rate.  The  bottleneck  server  is  assumed  to  be  highly 
sensitive  to  errors,  and  a  large  weight  is  attached  to  this 
server  in  the  iteration  process  (see  Chapter  3).  To  obtain 
an  M/G/l/N  solution,  the  general  service  time  distribution 
is  characterized  by  its  first  two  moments  and  is  represented 
by  Cox-type  servers  as  discussed  in  Section  2. 

In  general,  using  more  information  about  the 
service  time  distribution  should  increase  the  accuracy  of 
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the  solutions.  The  EPF  solution  technique  does  well  in 
most  of  the  test  cases  presented  in  [Shum  1976]. 

Shum  makes  an  observation  about  the  accuracy  of  EPF 
based  upon  the  value  of  the  CV  of  the  general  service  time 
distribution;  this  is  called  the  EPF  and  C  Hypothesis  (C 
stands  for  the  CV) .  According  to  this  hypothesis,  EPF  does 
well  (in  terms  of  accuracy)  for  a  CV  close  to  1  and  for  very 
large  values  of  the  CV  (CV>400).  Similar  observations 
for  the  Norton  reductions  have  been  made  by  Tripathi  and 
Levy  [1976].  We  discuss  error  behaviour  in  greater  detail 
in  Chapter  5;  we  show  that  this  error  behaviour  is  not  due 
to  any  particular  approximation  technique,  but  is  due  to  the 
structure  of  the  Cox-type  server  used  to  represent 
non-exponential  distributions. 


NUMERICAL 

Recently,  Marie  [1979]  characterized  each  server  as 
a  load  dependent  server,  with  load  dependent  arrivals. 
Thus,  each  centre  is  assumed  to  be  a  load  dependent  server 
with  general  service  time  distribution  and  with  load 
dependent  Poisson  arrivals. 

Since  absolute  arrival  rates  at  the  servers  are  not 
known,  an  iterative  technique  is  suggested.  The  iteration 
is  performed  until  desired  levels  of  accuracy  in  two 
invariants  are  achieved.  The  two  invariants  used  are:  the 

total  of  the  queue  lengths  at  all  the  centres  should  be 
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equal  to  the  number  of  customers 
throughput  constraint  should  be 
individual  centres  is  obtained 
Marie  has  reported  a  very  fas 
accuracy  for  this  approximation. 


in  the  network,  and  the 
satisfied.  The  solution  to 
by  numerical  techniques, 
t  convergence  and  very  good 


2.4.3  Other  Approximations 

In  this  section  we  discuss  two  types  of 
approximations:  approximations  based  on  mean  value  analysis, 
and  approximations  that  change  the  model  structure  in  order 
to  obtain  an  efficient  solution. 

Recently,  based  on  mean  value  analysis  [Reiser  and 
Lavenberg  1978,  Reiser  1979],  Bard  [1979]  has  proposed  an 
approximation  technique  to  solve  queueing  network  models 
that  do  not  possess  product  form.  Following  Bard's 
notation,  let  Nij,  Tij,  Fij,  and  N(ij)  be  the  average  number 
of  class  i  customers  in  queue  j,  the  average  time  spent  by 
class  i  customers  in  queue  j,  a  function  whose  form  depends 
on  the  queueing  discipline,  and  the  average  conditional 
queue  population  at  the  time  that  a  class  i  customer 
arrives  at  queue  j,  respectively.  The  network  has  M  centres 
and  K  classes.  Also,  let  ni  be  class  i  population  size  and 
dij  be  class  i  throughput  at  queue  j.  Then  by  Little's  law 

Nij  =  dij*Tij  (2.4) 


Also , 
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ni  =  Nil+Ni2+. . .+NiM  (2-5) 

and 

Tij  =  Fij (N ( i j ) )  (2.6) 

Summing  (2-4)  over  j  and  using  (2-5),  we  get 

Nij  =ni*Ti j/(Til+Ti2+. - -+TiM)  (2-7) 

Bard  assumes  that  N(ij)  is  the  same  as  the  distribution  of 
Nij  in  the  steady  state  and  proposes  the  following  iterative 
algorithm  to  solve  for  the  steady  state: 

1-  Assume  an  initial  value  for  Nij. 

2-  Compute  Tij  using  (2-6). 

3-  Recompute  Nij  using  (2-7). 

4-  Return  to  step  2- 

The  iteration  is  continued  until  no  significant  change  in 
Nij  from  one  iteration  to  the  next  is  observed. 

Reiser  [1976]  proposed  an  approximation  technique 
for  solving  queueing  network  models  with  priority  scheduling 
disciplines-  He  treats  a  model  involving  two  priority 
classes.  He  suggests  a  decomposition  of  the  model  in  which 
the  lower  priority  class  is  effectively  served  by  a 
processor  whose  capacity  is  reduced  by  the  utilization 
attained  by  high  priority  jobs. 

Sevcik  [1977]  provided  upper  and  lower  bounds  for 
the  performance  measures  of  a  similar  two  class  system  with 
priority.  The  major  differnce  is  that  in  Sevcik's  model, 
the  two  classes  contend  with  one  another  at  the  I/O  devices 
(on  a  FCFS  basis)  as  well  as  at  the  CPU  (on  a  priority 
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basis).  Sevcik  suggests  modifications  to  the  model  that 
make  efficient  exact  solution  possible  and  that  yield 
answers  guaranteed  to  be  upper  or  lower  bounds  on  the 
performance  of  the  original  model. 

Such  approximations  as  these  change  the  model  into 
one  for  which  the  exact  solution  can  be  efficiently 
attained.  Because  the  methods  discussed  in  this  thesis  deal 
with  approximate  solutions,  they  are  not  relevant  to 
approaches  that  force  models  into  forms  for  which  efficient 
exact  solution  is  possible. 


A  UNIFIED  APPROACH  TO  IMPROVING  APPROXIMATIONS 


This  chapter  examines  ways  of  improving  the  approximate 
solution  techniques  discussed  in  Chapter  2.  We  show  that 
various  approximation  techniques,  although  different  in 
origin,  presentation,  and  motivation,  have  common  sources  of 
errors-  A  unified  approach  that  can  be  useful  to  improving 
the  accuracy  of  some  of  these  approximations  is  presented. 

Section  3-1  identifies  various  sources  of  errors  in 
both  aggregation/decomposition  and  product  form 
approximations.  Section  3.2  establishes  a  relationship 
between  the  sources  of  errors  in  both  types  of  approximation 
and  presents  a  unifying  hypothesis.  Section  3.3  uses  this 
unifying  hypothesis  to  show  that  a  single  effort  may  improve 
more  than  one  approximation.  In  particular,  we  show  that  a 
better  estimation  of  the  interactions  between  the  service 
centres  can  improve  both  the  Norton  theorem  reduction  and 
the  diffusion  approximation.  Examples  of  both  open  and 
closed  systems  are  presented.  Section  3.4  summarizes  the 
results  of  this  chapter. 
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3-1  Problems  In  Approximations 

To  improve  an  approximate  solution  technique,  we 
must  first  identify  why  the  approximations  are  introduced. 
To  do  this  we  must  examine  the  assumptions  made  in  arriving 
at  a  solution  to  the  model. 

For  example,  consider  the  interactive  system 
presented  in  Figure  2.6.  The  decomposition  technique  gives 
a  good  approximate  solution  if  there  is  only  limited 
interaction  between  the  two  subsystems  (the  terminals  and 
the  central  subsystem).  Therefore,  if  accurate  solutions  to 
the  individual  subsystems  can  be  obtained,  the  degree  of 
interaction  among  the  subsystems  largely  determines  the 
error . 

The  next  two  subsections  identify  various  sources 
of  errors  for  both  aggregation/decomposition  and  product 
form  approximations.  Section  3.1.1  examines  the  accuracy  of 
aggregation/decomposition  techniques.  In  Section  3.1.2 
product  form  approximations  are  discussed. 


3.1.1  Improving  Aggregation/Decomposition  Approximations 

Aggregation/decomposition  is  based  on  modeling  the 
behaviour  of  an  entire  subsystem  by  a  single  load  dependent 
composite  server.  The  service  time  distribution  of  the 
composite  server  is  specified  as  the  load  dependent 
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inter-departure  times  from  the  subsystem-  If  the  conditions 
for  a  product  form  solution  are  met  or  the  system  is 
completely  decomposable,  the  set  of  load  dependent  mean 
service  rates  characterizes  the  distribution  completely.  In 
this  case  the  composite  server  is  termed  a  true  flow 
equivalent  server  [Chandy  and  Sauer  1978]. 

When  the  conditions  for  a  product  form  solution  are 
not  met,  however,  specifying  a  true  flow  equivalent  server 
is  an  open  problem.  At  present,  the  best  we  can  do  is 
derive  an  approximate  flow  equivalent  server. 

The  two  principal  factors  that  characterize  an 
approximate  flow  equivalent  server  are:  (1)  the 
inter-departure  times,  and  (2)  effects  from  the  outside  of 
the  subsystem. 


INTER-DEPARTURE  TIME  DISTRIBUTION 


The  inter-departure  time  distribution  is  a  local 
characteristic  of  a  subsystem  that  depends  only  on  the 
number  of  customers  in  the  subsystem.  In 
nearly-decomposable  systems,  the  inter-departure  time 
distribution  is  asymptotically  exponential  as  the 
interactions  among  the  subsystems  approach  zero  [Chandy  and 
Sauer  1978].  Thus,  only  the  mean  inter-departure  times  are 
needed.  The  situation  is  different,  however,  in  Norton's 
theorem  aggregation.  It  is  this  case  that  we  address  in  the 
remainder  of  this  subsection. 
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Calculating  the  inter-departure  time  distribution 
is  in  general  a  hard  problem.  Moreover,  the  distribution 
depends  on  the  number  of  customers  in  the  subsystem.  Thus, 
even  if  we  had  the  exact  inter-departure  time  distribution 
for  each  possible  number  of  customers,  existing  methods 
cannot  properly  represent  this  in  the  overall  model.  (Note 
that  we  can  represent  load  dependent  mean  service  times,  but 
not  load  dependent  service  time  distributions.) 

Even  if  we  could  be  content  with  knowing  the  first 
two  moments  of  the  inter -departure  time  distribution  of  the 
aggregated  server,  problems  would  still  remain.  The  second 
moment  of  the  inter-departure  time  distribution  is  also 
dependent  on  the  load,  and  we  cannot  represent  such 
distributions  in  our  analysis.  Thus  having  a  method  to 
estimate  the  load  dependent  first  two  moments  of  the 
inter-departure  time  distribution  of  the  aggregated  server 
is  of  no  use. 

Calculating  the  load  dependent  mean  inter-departure 
time  involves  solving  the  subsystem  for  throughput  at 
various  loads.  The  most  accurate  solution  technique,  the 
global  balance  solution  method,  however,  is  also  the  most 
expensive  technique. 

One  way  to  reduce  the  complexity  of  solving  the 
subsystem  by  the  global  balance  technique  is  by  repeated 
application  of  aggregation  [Zahorjan  1977].  For  example, 
suppose  we  had  a  program  to  determine  efficiently  the  steady 
state  solution  of  a  network  with  j  service  centres. 
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Grouping  the  subsystem  into  further  subsystems  with  at  most 
j  service  centres,  each  subsystem  can  be  solved  separately 
by  the  global  balance  technique  and  is  represented  by  an 
aggregated  server-  This  aggregation  could  be  applied 
repeatedly  until  a  solution  to  the  subsystem  is  obtained. 
At  the  second  and  subsequent  stages  of  the  aggregation 
process,  however,  only  the  mean  inter-departure  time  of  the 
aggregate  server  is  known-  Thus,  at  each  stage  the  second 
moment  of  the  inter-departure  time  distribution  has  to  be 
approximated  -  propagating  the  error  in  estimating  the 
second  moment  at  each  subsequent  stage. 

An  alternative  approach  to  using  the  global  balance 
solution  technique  for  solving  the  inter-departure  time 
distribution  of  the  aggregated  server  is  to  assume  that  the 
subsystem  satisfies  the  conditions  of  product  form  and  solve 
it  by  the  product  form  solution  technique  [Sauer  and  Chandy 

i 

1975]-  This  method,  although  efficient,  can  introduce  large 
error  [Tripathi  and  Levy  1976]. 

To  approximate  the  load  independent  CV  of  the 
inter-departure  time  distribution  of  the  aggregate  server, 
equation  (2.1)  can  be  used.  There  are  some  basic  drawbacks 
in  this  simple  approximation  [Sevcik  et  al.  1977].  Further 
improvements  in  estimating  the  CV  are  discussed  in  Section 


3.  3. 
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OUTSIDE  EFFECTS 

The  principal  advantage  as  well  as  the  principal 
disadvantage  of  aggregation  is  that  each  subsystem  is  solved 
separately.  This  makes  the  problem  more  tractable,  but 
ignores  any  effects  due  to  interaction  among  subsystems. 
These  interactions  represent  the  global  behaviour  or  the 
outside  effect  on  the  subsystem. 

In  calculating  the  inter-departure  time 
distribution  for  Norton's  theorem  aggregation,  the  subsystem 
is  analyzed  off-line  (i.e.,  the  subsystem  assumes  no 
knowledge  of  the  environment  in  which  it  is  supposed  to 
work).  This  method  assumes  that  the  output  process  from  the 
subsystem  is  the  same  as  the  input  process  to  the  subsystem. 
In  general,  however,  these  two  processes  are  different. 
Therefore,  assuming  that  they  are  identical  can  introduce 
error.  A  way  of  reducing  this  error  is  presented  in  Section 
3.3. 


3.1.2  Improving  Product  Form  Approximations 


If  a  system  has  a  product  form  solution,  the 
solutions  for  the  service  centres  are  independent. 
Moreover,  each  service  centre  is  completely  characterized  by 
the  mean  arrival  rate  and  the  mean  service  time.  In  open 
systems,  mean  arrival  rates  to  service  centres  can  be 
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calculated  explicitly-  In  closed  systems,  these  arrival 
rates  can  only  be  calculated  up  to  a  multiplicative 
constant,  but  absolute  arrival  rates  are  not  necessary  for 
product  form  systems. 

For  the  systems  that  do  not  satisfy  product  form 
conditions,  however,  the  arrival  process  has  to  be 
calculated  explicitly-  Using  relative  arrival  rates  can 
lead  to  error  in  the  solution  to  the  network  [Shum  1976] - 
Thus,  in  characterizing  errors  in  product  form 
approximations,  two  factors  have  to  be  considered:  (1) 
obtaining  solutions  to  individual  service  centres,  and  (2) 
estimating  the  absolute  arrival  rates  at  individual  servers. 


SOLUTION  TO  INDIVIDUAL  SERVICE  CENTRES 

The  principal  approaches  involving  product  form 
approximations  are  EPF  [Shum  1976]  and  diffusion  [Kobayashi 
1974].  EPF  assumes  that  the  arrival  process  at  each  centre 
is  Poisson.  Each  centre  is  considered  as  an  M/G/l/N  system. 
Solutions  to  M/G/l/N  systems  are  derived  analytically 
assuming  the  knowledge  of  only  the  first  two  moments  of  the 
service  time  distribution  and  the  distributional 
representations  discussed  in  Chapter  2.  Thus,  there  is  no 
local  error  introduced  in  EPF.  However,  error  is  introduced 
by  the  assumption  that  the  service  time  distribution  is 
assumed  to  have  a  specific  form.  We  view  this  error  as 
modeling  error.  It  is  important  to  distinguish  such  error 
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arising  from  modeling  assumptions  from  the  error  due  to 
inexact  solution  of  a  model. 

In  diffusion  approximations,  solving  each 
individual  centre  involves  many  approximations.  The  service 
time  distribution  is  characterized  only  by  the  first  two 
moments.  The  queue  length  is  represented  by  a 
continuous-valued  variable  which  requires  mapping  continuous 
queue  lengths  into  discrete  integral  queue  lengths.  The 
errors  introduced  by  these  assumptions  have  already  been 
discussed  in  Chapter  2. 

ESTIMATING  THE  ARRIVAL  PROCESSES 

EPF  approximations  assume  that  the  arrivals  at  all 
the  servers  are  Poisson.  Thus,  only  the  mean  arrival  rate 
is  needed.  In  open  systems,  the  mean  arrival  rate  is  easily 
obtainable,  but  not  in  closed  systems. 

Absolute  arrival  rates  are  calculated  by  an 

iterative  process  that  uses  the  fact  that  if  Xl,X2,«..,XM 

are  the  absolute  arrival  rates  then  they  must  satisfy  the 

throughput  constraint,  namely 

M 

Xi  =  ^  \  j  *  P  j i  For  i  =  l,2,.«.,M 

where  Pji  is  the  transition  probability  from  centre  j  to  i. 
Assuming  XI, X2, . . . . , XM  as  absolute  arrival  rates,  we 
calculate  server  throughputs  using  the  EPF  solution 
technique.  Let  Tl, T2, . . . . , TM  be  the  resulting  throughputs. 
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Then  Tl, T2, - - . - , TM  should  also  satisfy  the  throughput 
constraint-  Therefore, 

(Xl/Tl)  =  (X2/T2)  = _ =  (XM/TM) 


Since  a  product 

form 

solution  is 

an 

approximation,  this 

equality  in  general 

may  not 

hold . 

Dur 

ing  each  iteration,  a 

new  arrival  rate 

at 

each 

centre 

is 

calculated  and  the 

iteration  is 

continued 

until 

the 

above  equation  is 

approximately  satisfied-  The  solution  thus  obtained  may  not 
be  the  exact  arrival  rates.  Although  the  convergence  of 
this  iterative  technique  is  not  guaranteed,  Shum  has 
reported  good  convergence  for  the  examples  examined-  For 
some  cases,  however,  there  may  be  some  problems  in 
convergence  [Balbo  and  Denning  1979]  - 

For  diffusion  approximations  the  arrival  process  at 
a  centre  is  approximated  using  estimates  of  the  first  two 
moments  of  inter-arr ival  times-  To  estimate  the  CV  of  the 
arrival  process  the  mean  arrival  rate  is  needed  first. 

Two  approaches  to  estimating  the  mean  arrival  rate 
have  been  suggested:  assume  that  the  bottleneck  server  is 
saturated,  or  solve  the  system  with  the  assumption  that 
product  form  exists-  Both  of  these  methods  are  quick  and 
easy  ways  to  estimate  the  mean  arrival  rates. 

Three  approaches  have  already  been  discussed  in 
Section  2.4.2  for  estimating  the  CV  of  the  arrival 
processes.  The  accuracy  of  these  estimates  depends  upon  the 
accuracy  of  the  mean  arrival  rates  and  upon  the  accuracy  of 
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the  formulas  both  for  the  splitting  and  merging  of 
non-Poisson  processes  and  for  the  departure  process  of  a 
G/G/l  system. 

In  short,  EPF  puts  a  great  deal  of  emphasis  on 
approximating  the  mean  arrival  rate  as  accurately  as 
possible.  Diffusion  approximations,  on  the  other  hand,  get 
quick  and  possibly  inaccurate  estimates  for  the  mean  arrival 
rates  and  devote  a  great  deal  of  effort  to  the  CV's  of  the 
arrival  processes. 


Having  examined  the  sources 

of 

err 

or  s 

in  various 

approx 

imate  solution  techniques, 

we 

are 

now 

ready 

to 

es tabl 

ish  a  close  relationship  between 

the 

two 

kinds 

of 

approximations. 

3.2  Unified  Approach  To  Solutions 

The  various  approximate  solution  techniques  can  be 
viewed  as  different  ways  of  achieving  a  simple  goal: 
breaking  a  complex  system  into  simpler  subsystems  that  can 
be  solved  efficiently  and  then  combining  these  subsystem 
solutions  to  get  the  overall  solution.  This  unified  view 
allows  the  errors  introduced  in  such  techniques  to  be 
divided  into  two  classes:  global  and  local. 

In  solving  a  subsystem,  global  error  is  introduced 
when  the  behaviour  of  the  rest  of  the  system  is  approximated 
(decomposition)  ;  local  error  is  introduced  when  the 
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subsystem  itself  is  approximated  (aggregation). 

Both  aggregation/decomposition  and  product  form 
errors  fit  very  nicely  into  the  above  classification.  The 
global  error  in  decomposition  is  due  to  interactions  between 
the  subsystems.  In  Norton's  theorem  approximations,  global 
error  results  from  the  fact  that  short  circuiting  a  part  of 
the  system  simply  ignores  global  interactions.  The 
decomposition  method  assumes  that  exact  solutions  to  the 
subsystems  can  be  obtained  and  thus  that  there  is  no  local 
error.  Norton's  theorem  approximations,  on  the  other  hand, 
can  introduce  local  error  -  either  by  assuming  product  form 
when  it  does  not  hold  or  by  repeated  applications  of 
aggregation . 

For  product  form  solutions,  each  centre  can  be 
considered  as  a  subsystem  we  are  trying  to  solve  separately. 
Global  error  is  introduced  in  EPF  approximation  due  to  the 
assumption  that  the  arrivals  at  each  subsystem  (centre)  are 
Poisson,  and  in  diffusion  approximations  due  to 
approximations  to  the  general  arrival  process.  Local  error 
is  introduced  in  EPF  and  diffusion  approximations  because 
only  the  first  two  moments  are  used  to  represent  the  general 
service  time  distributions.  Further  error  is  introduced  in 
diffusion  approximations  because  G/G/l  solutions  are 
approx imate . 


In  Bard's  approximation  [1979],  based  on  mean  value 
analysis,  global  error  is  introduced  by  assuming  that  N(ij), 
the  distribution  of  queue  population  at  the  time  that  a 
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class  i  customer  arrives  at  queue  j,  is  the  same  as  the 
steady  state  distribution  of  the  number  of  customers,  when 
calculating  Tij,  the  average  time  spent  by  class  i  customers 
in  queue  j-  Local  error  is  incurred  due  to  the  fact  that 
only  the  mean  Tij  is  considered  in  the  computations-  Thus, 
arrival  processes  are  characterized  by  only  their  mean. 

This  classification  has  two  basic  advantages  in 
studying  the  accuracy  of  approximate  solution  techniques: 

1-  It  increases  our  understanding  of  errors  in 

approximations  by  identifying  and  relating  them  to  the 
problems  in  other  approximations.  For  example,  a  good  way 
of  removing  global  error  in  one  approximation  technique  can 
also  be  used  to  improve  global  error  in  a  different 

technique.  An  illustration  of  this  is  given  in  the  next 
section. 

2.  It  suggests  new  approximation  techniques  by  identifying 
sources  of  error.  One  such  new  technique  is  presented  in 
Chapter  4. 

The  various  approximations  can  be  characterized 
based  upon  the  following  four  considerations: 

(a)  service  time  distribution  at  centres, 

(b)  transition  (inter-event)  processes, 

(c)  individual  subsystems  and  their  solutions,  and 

(d)  combining  subsystem  solutions  to  obtain  the  overall 
solution . 

The  characterizations  are  presented  in  Table  3.1. 
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Table  3.1 

1.  Product  Form  [Baskett  et  al.  1975] 

(a)  Phase  type 

(b)  Relative  arrival  rates 

(c)  Single  centres,  exact  Markovian  analysis 

(d)  Product  form  with  normalization  constant 

2.  Norton's  Theorem  Approximations 

(i)  By  Sauer  and  Chandy  [1975] 

(a)  Mean  and  CV 

(b)  Load  dependent  Poisson 

(c)  Arbitrary  subsystems,  product  form  approximations 

(d)  Global  balance  solution  with  flow  equivalent 

servers 

(ii)  By  SLTZ  [Sevcik  et  al.  1977] 

(a)  Mean  and  CV 

(b)  Mean  and  CV 

(c)  Arbitrary  subsystems,  product  form  approximations 

(d)  Global  balance  solution  with  flow  equivalent 

servers 

3.  Chandy  et  al.'s  [1975b]  Iterative  Approximation 

(a)  Mean 

(b)  Load  dependent  Poisson 

(c)  Single  centres,  M/G/l/N  solutions 

(d)  Iterative  approach  to  satisfy  invariants 

4.  EPF  [Shum  and  Buzen  1977] 

(a)  Mean  and  CV 

(b)  Poisson 

(c)  Single  centre,  M/G/l/N  analytic  solutions 

(d)  Product  form  with  normalization  constant  with  an 
iterative  technique  to  satisfy  throughput  constraint 

5.  Diffusion  Approximations  [Reiser  and  Kobayashi  1974] 

(a)  Mean  and  CV 

(b)  Mean  and  CV 

(c)  Single  centre,  G/G/l  diffusion  approximation 

(d)  Product  form  with  normalization  constant 

6.  Kuhn's  Approximation  [1976] 

(a)  Mean  and  CV 

(b)  Mean  and  CV 

(c)  Single  centre,  G/G/l  numerical  approximation 

(d)  Nothing  (only  open  systems  are  considered) 

7.  Marie's  [1979]  Approximation 

(a)  Phase  type 

(b)  Load  dependent  Poisson 

(c)  Arbitrary  subsystems,  numerical  approximation 

(d)  Iterative  procedure  to  satisfy  invariants 
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3-3  Examples  Of  Unified  Improvements 

This  section  demonstrates  the  benefits  of  the 
unified  approach  taken  in  this  chapter.  We  show  that  some 
of  the  techniques  developed  for  the  improvement  of  the 
aggregation/decomposition  approximations  can  also  be  used  to 
improve  product  form  approximations.  In  particular,  the 
characterization  of  the  global  behaviour  of  subsystems  in 
Norton's  theorem  aggregation  can  be'  used  to  reduce  the 
global  error  of  the  individual  service  centres  in  diffusion 
approximations.  The  remainder  of  this  section  presents  the 
improvement  methodology  and  gives  examples  to  demonstrate 
the  usefulness  of  this  methodology. 


3.3.1  Reducing  Global  Error 

After  presenting  some  improvements  for  Norton's 
theorem  aggregation,  this  section  discusses  how  these 
improvements  can  be  applied  to  diffusion  approximations. 
The  first  part  of  the  presentation  is  based  upon  results 
reported  in  [Sevcik  et  al.  1977]. 

IMPROVING  NORTON'S  THEOREM  AGGREGATION 

One  way  to  reflect  the  outside  interference  (global 
behaviour)  of  a  subsystem  in  Norton's  theorem  aggregation  is 
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to  adjust  the  second  moment  of  the  inter -depar ture  time 
distribution  of  the  composite  server.  This  can  be 
accomplished  by  estimating  the  second  moment  of  the  service 
time  distribution  of  the  composite  server  such  that  for  a 
given  input  process  the  composite  server  produces  the  same 
output  process  as  does  the  subsystem.  For  example,  consider 
the  model  in  Figure  3.1: 


Rest  of  the 

A 

The 

- ? 

System 

*  } 

Subsystem 

B 

- 4 — - — 

Figure  3.1:  Aggregation 


Given  a  transition  process  at  A,  the  composite 
server  should  yield  the  same  transition  process  at  B  as  the 
subsystem  does.  In  particular,  estimates  for  the 

transition  processes  at  A  and  B  are  obtained  and  the  CV  of 
the  composite  server  is  chosen  such  that  transition  process 
A  is  transformed  by  the  subsystem  into  approximately 
transition  process  B. 

Before  we  can  obtain  estimates  for  the  transition 
processes  at  A  and  B,  there  are  three  problems  to  address: 
the  splitting  of  a  transition  process,  the  merging  of 
transition  processes,  and  the  transformation  of  an 
transition  process  by  passing  through  a  server. 
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We  characterize  a  transition  process  by  the  mean 
and  the  variance  (or,  equivalently,  the  CV)  of  the 
inter-event  times.  The  transition  processes  are  assumed  to 
be  independent  of  each  other.  In  general,  the  transition 
processes  are  not  renewal;  however,  we  assume  that  they  are 
renewal  to  keep  the  problem  tractable.  Also,  successive 
inter-event  times  are  assumed  to  be  independent.  It  must  be 
noted  here  that  characterizing  a  transition  process  only  by 
the  first  two  moments,  assuming  independence  among  the 
processes,  and  assuming  them  to  be  renewal,  all  contribute 
error  in  the  analysis. 


SPLITTING  AND  MERGING  OF  PROCESSES 

The  splitting  and  merging  of  Poisson  streams  do  not 
present  any  problem;  both  these  operations  preserve  the 
Poisson  property.  The  situation,  however,  is  different  for 
non-Poisson  streams.  First,  results  are  derived  for 
splitting  a  stream  into  two  streams  and  merging  two  streams; 
the  results  for  the  general  case  are  obtained  later  by  the 
repeated  applications  of  the  results  for  two  streams.  The 
order  in  which  the  pairings  are  done  has  no  effect  on  the 
results  in  both  cases  (see  Appendix  A  for  a  discussion). 

Let  us  consider  splitting  a  stream  with  the  mean 
and  the  CV  of  the  inter-event  times  equal  to  Mx  and  CVx 
respectively.  Let  the  stream  be  split  by  sending  each  event 
one  of  the  two  ways  according  to  a  Bernoulli  trial  with 
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respective  probabilities  PI  and  P2  (Pl+P2=l),  and  Myl  and 
My 2,  and  CVyl  and  CVy2  be  the  means  and  the  CV's  of  the 
resulting  streams.  (See  Figure  3.2) 


Figure  3.2:  Splitting  a  Stream 


It  can  be  shown  that  (see  Appendix  A) : 

Myi  =  Mx/Pi,  and 

CVyi  =  1  +  Pi* (CVx-1)  for  i=l,  2  (3.1) 

It  can  be  seen  from  (3.1)  that  the  CV's  of  the 
resulting  processes  are  never  further  from  one  than  the  CV 
of  the  original  process  was.  Also,  the  CV  of  the  output 
process  selected  with  probability  P  asymptotically 
approaches  one  as  P  approaches  zero. 

The  first  two  moments  of  the  original  process 
completely  determine  the  first  two  moments  of  the  resulting 
processes.  In  fact,  it  can  be  proven  that  for  any  fixed 

n,  the  first  n  moments  of  the  original  process  are 
sufficient  to  calculate  the  first  n  moments  of  the  resulting 
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processes  (see  Appendix  A  for  the  proof) . 

Let  us  now  consider  the  merging  of  two  streams  with 
means  Mxl  and  Mx 2  and  CV’s  CVxl  and  CVx2  respectively.  Let 
the  resulting  process  be  represented  by  the  mean  My  and  the 
coefficient  of  variation  CVy.  (See  Figure  3.3) 


(Mxl,  CVxl) 


(Mx2, 


*(My,  CVy) 


The  merging  of  general  renewal  processes,  even  with 
no  serial  correlation,  is  more  complex.  The  expression  for 
My  can  be  obtained  by  simple  manipulation: 

My  =  (Mxl*Mx2)/ (Mxl+Mx2) 

The  closed  form  expression  for  CVy,  however,  is  not 
known,  and  complicated  integrals  are  involved  [Sevcik  et  al. 


1977,  Kuhn 

1976] . 

Kobayashi  [1974]  proposed 

an  approximate 

expression 

for 

CVy  (see  (2.2)).  In  the 

present  context 

(2.2)  reduc 

es  to 

n 

CVy  = 

1  X 

*(My/Mx i )  *  (CVxi  -  1) 

(3.2) 

3 

L-l 
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There  are  some  basic  shortcomings  in  using  (3.2). 
For  example,  consider  a  situation  when  CVxi  is  the  same  for 
all  i.  Then  CVy=CVxi.  It  is  known,  however,  that 
merging  non-Poisson  streams  yields  a  process  much  closer  to 
Poisson  than  the  original  ones.  The  following  approximation 
is  intended  to  remedy  this  situation. 


CVy 


1  +  ^  J My /Mx if* (CVxi 


1) 


(3.3) 


In  the  absence  of  analytic  results  about  CVy,  we 
compare  approximations  (3.2)  and  (3.3)  with  simulation 
results.  Two  examples  are  presented.  Example  1 

assumes  that  there  are  n  transition  processes  such  that 
Mxi=Mxj,  for  i , j =1 , 2, . . . , n ,  and  CVxi=9.0.  For  varying  n, 

the  results  obtained  by  (3.2),  (3.3),  and  simulations  are 
presented  in  Figure  3.4.  As  expected,  (3.3)  shows  a  very 
good  agreement  with  simulation  results.  Approximation 
(3.2),  however,  shows  a  large  error  with  respect  to 

simulation  results. 

In  the  second  example,  we  consider  two  processes 
with  CVxl=9.0  and  CVx2=2.0.  The  weights  of  xl  and  x2  are 
changed,  i.e.,  (My/Mxl)  is  varied  from  0  to  1 

(0, 1/4, 1/3, 1/2, 2/3, 3/4, 1) .  The  results  are  presented  in 
Figure  3.5.  Here  again,  (3.3)  has  better  accuracy  than 


(3.2). 


The  two  sets  of  examples  presented  here  and  the 
facts  that  (1)  merging  two  streams  with  CV  less  than  one 
will  a  yield  a  stream  with  somewhat  larger  CV,  (2)  merging 
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streams  with  CV  more  than  one  will  yield  a  stream  with 
somewhat  smaller  CV,  and  (3)  merging  arbitrary  large  numbers 
of  streams  with  the  same  CV  will  yield  a  nearly  Poisson 
stream,  all  support  the  contention  that  (3.3)  is  usually 
better  than  (3.2). 

For  splitting  and  merging  of  more  than  two  streams 
both  (3.2)  and  (3.3)  generalize  straightforwardly  (Appendix 
A).  In  the  rest  of  the  thesis  equations  (3.2)  and  (3.3)  will 
be  referred  to  in  their  generalized  terms,  treating  more 
than  two  streams. 

In  closing  we  must  again  emphasize  the  fact  that  in 
all  the  above  discussion  it  is  assumed  that  the  transition 
processes  are  independent.  This  is  hardly  the  case  in 
queueing  network  models.  Therefore,  the  results  obtained 
here  should  be  used  with  caution. 

PASSAGE  THROUGH  A  SERVER 

Given  a  general  arrival  process  (Ma,CVa)  and  a 
service  centre  with  general  service  time  distribution 
(Ms , CVs ) ,  we  want  to  determine  the  departure  process 
(Md,CVd).  In  the  steady  state  Md=Ma«  Therefore,  we  are 
interested  in  calculating  CVd.  (See  Figure  3.6) 
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(Ma.  CVa) 

(Ms,  CVs) 

— — - — > 

.(Md,  CVd) 


Figure  3.6:  Passage  Through  a  Server 


For  M/G/l  systems  (that  is,  CVa=l)  it  has  been 
proven  by  Disney  and  Cherry  [1974]  that,  if  u=Ms/Ma  then 

CVd  =  1  +  ui*(CVs  -  1) 


There  is  no  such  closed  form  expre 
arrival  process  is  non-Poisson, 
process  depends  heavily  upon  the 
server  is  lightly  loaded  (i.e. 
heavily  on  the  service  process  if 
loaded . 


ssion  for  CVd  when  the 
Observe  that  the  departure 
arrival  process  if  the 
,  u  is  small)  and  depends 
the  server  is  heavily 


Gelenbe  and  Pujolle  [1976]  suggested  (2.3)  as  an 
approximation  to  the  departure  process  of  a  G/G/l  server. 
Their  method  produces  incorrect  answers  for  some  cases  where 
exact  results  are  known.  For  example,  (2.3)  yields  a  value 
of  u(l-u)  for  the  D/D/1  case,  whereas  the  exact  value  is 
zero  [Pack  1975]. 


Results  based  on  simulations  for  a  range  of  values 
of  u,  CVa ,  and  CVs  along  with  the  observation  about  the 
dependence  of  the  departure  process  on  the  arrival  and 
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service  processes  lead  to  the  following  approximation: 

CVd  =  1  +  (1-u2)  *  (CVa-1)  +  ua*  (CVs-1)  (3.4) 

Kuhn  [1976]  independently  reported  an  expression  identical 
to  (3.4);  his  result  is  based  upon  approximating  the  G/G/l 
expression . 

The  expression  (3.4)  is  exact  for  M/G/l  and  D/D/1 
systems.  Figure  3.7  presents  a  comparision  between  the 
results  obtained  by  the  simulations,  equation  (2.3),  and 
equation  (3.4).  The  simulation  program  uses  the 
distributional  forms  dicussea  in  Figure  2.3  to  represent 
various  values  of  CVa  and  CVs.  Equation  (3.4)'s  predictions 
are  significantly  closer  to  the  simulation  results  than 

(2.3) 's«  Additionally,  Table  3.2  compares  results  obtained 
by  equations  (2.3)  and  (3.4)  in  four  specific  cases.  For 
the  M/M/1  and  M/D/1  cases  both  of  them  agree  and  give  exact 
results  [Muntz  1973,  Pack  1975].  For  the  D/M/1  case  (2.3) 
and  (3.4)  give  different  results  although  they  approach  one 
another  as  the  utilization  approaches  zero  or  one.  For  the 
intermediate  utilizations,  simulation  results  verify  that 

(3.4)  gives  more  accurate  estimate  of  CVd  than  (2.3). 
Finally,  for  the  D/D/1  case  equation  (3.4)  gives  the  correct 
result  while  (2.3)  does  not. 
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Figure  3.7:  Departure  Processes 
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Table  3.2 


Case  CVa  CVs 

M/M/1  1  1 

M/D/1  1  0 

D/M/1  0  1 

D/D/1  0  0 


CVd 

Eq  -  (  2- 3 )  Eq «  (3.4) 
1  1 

l-(u*u)  1- (u*u ) 

u  (u*u) 

u* ( 1-u )  0 


To  combine  the  three  results  discussed  above  and  to 
demonstrate  their  utility,  let  us  consider  an  example  of  a 
closed  central  server  model  (Figure  2.1).  Suppose  that  we 
want  to  approximate  the  I/O  subsystem  by  an  aggregate 
server.  Estimates  for  the  mean  arrival  rates  at  the  servers 
are  needed  before  we  can  use  the  results  of  (3.1),  (3.3), 
and  (3.4).  We  calculate  a  first  approximation  of  the 
central  server  utilization  using  the  basic  Norton's  theorem 
decomposition  described  in  Chapter  2.  Having  obtained  the 
central  server  utilization,  the  mean  arrival  rates  at 
various  servers  can  be  calculated.  Let  Mai  be  the  mean 
inter-arrival  time  at  the  central  server.  Based  upon  the 
results  of  (3.1),  (3.3),  and  (3.4),  the  CV  of  the  arrival 
process  at  the  central  server  is  (see  Appendix  B  for  the 
der ivation ) : 
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(3.5) 


where  Msl  and  CVs 1  are  the  mean  and  the  CV  of  the  service 
time  distribution  at  the  central  server,  i=2,...,M  are  the 
I/O  servers  whose  service  time  distributions  have  means  Msi 
and  coefficients  of  variation  CVsi,  and  Pi  is  the  transition 
probability  from  the  central  server  to  the  i-th  server. 

Given  CVal,  CVdl  can  be  obtained  by  using  (3.4). 
Thus,  we  can  estimate  the  arrival  process  to  the  subsystem, 
(Mai, CVdl),  and  the  departure  process  from  the  subsystem, 
(Mai, CVal)  (recall  that  Mal=Mdl).  The  CV  of  the 
inter-departure  time  of  the  aggregate  server  is  adjusted  so 
that  the  server  yields  the  same  transition  processes.  Sevcik 
et  al.  [1977]  have  shown  that  this  technique  reduces  the 
error  in  the  performance  measures  obtained  by  the  Norton's 
theorem  approximation  in  some  cases.  We  will  refer  to  this 
approximation  as  Norton's  theorem  by  SLTZ  from  now  on. 

In  a  general  network,  the  CV's  of  the  transition 
processes  can  be  calculated  by  solving  a  set  of  linear 
equations  representing  the  relationships  among  the  departure 
processes  of  the  various  service  centres.  The  departure 
process  at  the  i-th  server  is  given  by  (see  Appendix  B) : 
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cvdi  -  1  ♦  fe)2*(cvsi  -  i,  ♦  d-(^)2)  E  fe)2^i3Mcvdj  -  D 

'  j  =  l  * 


(3.6) 


MODIFIED  DIFFUSION 

In  order  to  calculate  the  diffusion  solution  at 
each  server,  we  need  the  mean  and  the  CV  of  the  arrival 
process  at  each  server.  The  results  presented  for  the 
global  improvement  of  the  Norton's  theorem  aggregation  can 
be  used  in  estimating  these  parameters. 

Using  approximation  (3.2)  for  merging  event  streams 
and  approximation  (2.2)  for  departure  processes  of  G/G/l 
systems  one  can  get  an  equation  similar  to  (3.6).  (See 
Appendix  B) 

CVdl  ■  1  .  (g4)!-(CV,.  -  1)  .  11  - 

- »] 

Gelenbe  and  Pujolle's  [1976]  method  implies  (3.7) 
is  a  good  estimate  for  CVdi.  They  showed  that  they  could 
improve  over  the  original  approximation  of  Kobayashi.  Since 
(3.3)  is  better  than  (3.2)  for  estimating  the  CV  of  the 
merged  process  and  (3.4)  is  better  than  (2.2)  for  estimating 
the  output  process  of  a  general  server,  these  improvements 
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can  be  used  to  improve  on  existing  diffusion  approximations. 

Thus,  we  can  improve  the  diffusion  approximations 


using 

the 

techniques  developed 

for  the 

aggregation 

approximation 

by 

using 

equation  (3.6) 

instead 

of 

equation 

(3.7). 

In 

the 

next 

section,  improvements 

in 

diffusion 

approximations  are  shown  by  using  this  suggestion. 


3.3.2  Examples 


We  present  examples  of  queueing  network  models  of 
computer  systems  to  show  that  the  techniques  discussed  in 
Section  3.1  improve  previous  approximations.  First  we 
present  open  network  examples  and  then  closed  networks. 


OPEN  SYSTEMS 


The  model  shown  in  Figure  3.8  represents  an 
interactive  computer  system.  This  model  has  been  used  by 
Reiser  and  Kobayashi  [1975]  and  also  by  Gelenbe  and  Pujolle 
[1976]  to  demonstrate  the  effectiveness  of  their 
approximations.  We  shall  use  their  parameters  in  order  to 
compare  our  improvement  with  their  approach. 

Jobs  arrive  to  the  system  according  to  a  Poisson 
process  with  mean  inter-arrival  time  equal  to  MaO.  After 
passing  through  the  first  server,  jobs  leave  the  system  with 
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probability  1-Pl  or  enter  the  queue  for  the  second  server 
with '^probability  PI-  After  service  at  the  second  server,  a 
job  may  enter  the  queue  at  server  1  with  probability  P2  or 
rejoin  the  queue  at  server  2  with  probabilty  1-P2-  Servers 
1  and  2  can  be  used  to  represent  a  CPU  and  an  I/O  device, 
respectively. 

Let  PI  and  P2  both  be  equal  to  0.5  and  MaO  be  equal 
to  2.  We  compare  the  results  obtained  by  various 
approximations.  The  approximations  suggested  by  Kobayashi 
[1974],  by  Gelenbe  [1976],  and  modified  diffusion  (based  on 
equation  (3.5))  are  compared  with  the  results  obtained  by 
simulations  and  by  assuming  exponential  service  time 
distributions  at  both  the  servers.  The  results  are 
presented  in  Tables  3.3  and  3.4.  Errors  are  expressed 
relative  to  the  simulation  results,  which  are  inexact. 
Although  enough  transactions  were  simulated  to  yield  a  95% 
confidence  interval  of  width  less  than  one  percent  of  the 
estimated  value,  results  must  be  interpreted  with  care. 
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Table  3.3 

Ms  1  =  1/. 9 ,  CVs  =  . 5 ,  Ms  2  =  1/. 84,  CVs2  =  0 


Method 

Queue 

1 

Queue 

2 

Mean  Length 

Rel . 
Error 

Mean  Length 

Rel. 

Error 

Exponential 

9.00 

30% 

5.  25 

65% 

Kobayashi 

6.76 

1.  5% 

2.  70 

15% 

Gelenbe  and 

Pu  jolle 

6.65 

2.5% 

2.95 

8% 

Equation  (3.5) 

Simulation 

6.82 

6.84 

.  2% 

2.96 

3.  22 

8% 

Table  3. 

4 

Ms  1=1/. 95,  CVs  = 

.5,  Ms 2 

=1/. 89 ,  C Vs  2  = 

0 

Method 

Mean  Length 

Rel. 

Error 

Mean  Length 

Rel. 

Error 

Exponential 

19.00 

43% 

9.00 

100% 

Kobayashi 

14.  30 

7.  5% 

4.  52 

1% 

Gelenbe  and 

Pu jolle 

12.  62 

5% 

4.  39 

2.  5% 

Equation  (3.5) 

12.83 

3.  5% 

4.40 

2.  5% 

Simulation 

13.  30 

— - 

4.  50 

Tables  (3. 

3)  and  (3.4) 

show  that  a  better 

estimate 

the  arrival 

processes  at 

the 

servers  can 

improve 

diffusion  approximation. 
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CLOSED  SYSTEMS 

We  consider  the  central  server  model  shown  in 
Figure  3.9.  Let  the  mean  and  the  CV  of  the  service  time 
distribution  at  server  i  be  Mi  and  CVi.  Let  N  be  the  number 
of  jobs  present  in  the  system. 

For  small  values  of  N  (N<9)  QSOLVE  [Levy  1977]  is 
used  to  obtain  exact  solutions  using  the  global  balance 
solution  technique.  For  larger  N,  however,  this  solution 
package  becomes  very  expensive,  and  we  have  to  rely  on 
simulation  results  to  compare  the  accuracies  of  the 
approximations. 


Le 

t 

P2 

=P3  = 

.3, 

Tl 

n 

» 

Ml 

=  4, 

CVI 

=Varying 

M2 

n 

o 

CV2 

=  9 

M3 

H 

i-* 

NJ 

CV3 

=1 

M4 

II 

I-1 

O 

CV4 

=  16 

N  = 

Varying 

Two  sets  of  experiments  are  performed:  for  the  first  set  of 
experiments,  N=3  and  CVI  is  varied  (4,9,25,64,100,400);  for 
the  second  set  of  experiments  CV1=4  and  N  is  varied 
(3,4, 5,6,7). 

Figures  3.10  and  3.11  present  the  error  behaviour 
for  the  utilization  of  the  first  server  for  three 
approximation  methods.  It  is  interesting  to  note  that 
although  the  suggested  approximation  improves  over  the 
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Arrivals 


Figure  3.8 


Figure  3. 9 


7o  Relative 

Error  in  CPU  Utilization 


Figure  3.10:  Error  in  CPU  Utilization  of  varying  C V  of  Service  Time  Distribution 
at  the  CPU 
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Figure  3.11:  Error  in  CPU  Utiliztion  for  Varying  Multi -programming  Level  (N) 
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existing  methods  in  most  of  the  cases  presented  here,  all 
these  approximations  get  worse  for  large  values  of  the  CV  of 
the  service  time  distribution.  The  main  objective  of  these 
examples  is  not  to  show  how  accurate  a  diffusion 
approximation  is,  but  to  demonstrate  that  improving  the 
estimates  for  the  global  behaviour  does  improve  the  existing 
approximations,  and  this  has  been  accomplished.  Further 
improvements  to  diffusion  approximations  are  discussed  in 
Chapter  4. 


3.4  Summary  and  Comments 


Most  approximate  solution  techniques  can  be  viewed 
as  decomposition  techniques:  the  system  is  divided  into 
subsystems  consisting  of  single  or  multiple  service  centres, 
each  subsystem  is  solved  separately,  and  then  these  separate 
solutions  are  combined  to  get  the  overall  solution.  Two 
kinds  of  errors  introduced  by  the  decomposition,  local  and 
global,  were  discussed. 

The  unified  approach  suggested  in  this  chapter 
helps  one  to  understand  the  approximations.  This  approach 
also  helps  to  improve  the  approximate  solution  techniques. 
A  method,  which  is  known  to  improve  Norton's  theorem 
aggregation,  has  been  shown  to  improve  diffusion 
approximations  as  well.  Further  use  of  such  an  approach 
will  be  discussed  in  Chapter  4. 
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Chapter  4 

GENERALIZED  PRODUCT  FORM 

Based  upon  the  unified  approach  presented  in  the 
previous  chapter,  we  propose  a  new  solution  technique  in 
this  chapter.  The  proposed  approximation  technique 
introduces  two  levels  of  improvements  by  separating  the 
errors  in  solving  a  system  into  two  parts.  Thus,  the  best 
possible  solutions  available  at  each  level  can  be  used  to 
yield  a  better  approximate  solution  technique  than  most  of 
the  existing  approximations.  We  shall  use  this  chapter  as  a 
guideline  to  developing  new  approximations. 

The  chapter  is  divided  into  four  sections.  The  first 
section  presents  motivation  for  the  new  approximation 
technique.  Section  2  introduces  a  generalized  product  form 
solution.  Two  examples  of  such  an  approximation  are 
presented  in  Section  3.  The  first  approximation  is  based 
upon  diffusion  approximation  and  the  second  approximation 
uses  simulation.  Examples  are  discussed  to  examine  the 
accuracy  of  the  new  approximation  techniques.  Finally,  in 
Section  4,  the  results  obtained  in  this  chapter  are 
summarized. 
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4.1  Introduction  and  Motivation 


A  good  decomposition  approximation  minimizes  both  local 
and  global  errors.  In  other  words,  each  decomposed 
subsystem  is  solved  as  accurately  as  possible,  and 
interactions  between  subsystems  are  estimated  and  accounted 
for  precisely. 

Reducing  error  at  one  level  does  not  necessarily  leave 
the  error  at  the  other  level  unaffected.  In  general,  the 
two  errors  are  highly  inter-related.  For  example,  consider 
a  decomposition  where  each  server  is  solved  separately. 
Assume  the  arrival  process  at  each  server  to  be  Poisson. 
Thus  each  server  can  be  solved  locally  quite  accurately 
using  M/G/l  solutions  [Shum  1976].  If  we  estimate  the 
arrival  process  at  service  centres  more  accurately,  a 
non-Poisson  arrival  process  may  result;  solutions  to  G/G/l 
systems  are  needed.  At  present  no  exact  closed  form 
solution  to  a  G/G/l  system  exists  and  one  must  use 
approximations.  Therefore,  although  the  error  in 
characterizing  arrival  processes  at  service  centres 
(interactions  between  centres)  has  been  reduced,  the 
solutions  to  individual  subsystems  may  introduce  higher 
error  (the  local  error). 

Thus,  in  proposing  a  new  approximation  we  have  to 
decide  whether  improving  the  accuracy  at  one  level  of 
approximation  increases  any  error  at  the  other  level.  If 
not,  this  should  then  lead  to  a  better  approximation.  In 


-87- 


summary,  in  proposing  any  approximation,  not 
accuracy  at  individual  levels  important, 
overall  effect  has  to  be  considered. 


only  is 
but  also 


the 

the 


4.2  Generalized  Product  Form 

Consider  a  queueing  network  that  is  divided  into  M 
subsystems.  A  Generalized  Product  Form  (GPF)  solution  to 
the  system  is  defined  as 

P(S1,S2, ... ,SM)  =  ( 1/G (N) )*F1(S1)*...*FM(SM)  (4.1) 

where  P (S 1, S2, . . . , SM)  is  the  steady  state  probability  that 
the  system  is  in  a  state  such  that  the  ith  subsystem  is  in 
the  state  Si,  G(N)  is  the  normalization  constant,  and  Fi(Si) 
is  the  solution  to  the  ith  subsystem.  To  obtain  the 
solution  to  the  individual  subsystems,  the  arrival  process 
to  each  subsystem  is  needed.  If  a  system  yields  a  product 
form  solution,  in  (4.1)  we  assume  that  the  arrival  processes 
to  subsystems  are  the  same  as  departure  processes  from  them. 
This  assumption,  as  will  be  shown  in  proposition  4.1, 
guarantees  that  (4.1)  yields  exact  solution  when  it 
satisfies  the  conditions  for  product  form  (see  Chapter  2). 

When  the  system  does  not  satisfy  the  conditions  for 
product  form,  the  arrival  processes  at  various  subsystems 
are  obtained  by  considering  the  effect  from  the  entire 
system.  Ey  short  circuiting  the  subsystem,  its  effect  on 
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the  rest  of  the  system  is  ignored  and  thus  error  may  be 
introduced  [Sevcik  et  al.  1977].  Techniques  to  estimate 
arrival  and  departure  processes,  discussed  in  Chapter  3,  can 
be  used. 

Since  (4.1)  has  product  form  format,  the  convolution 
algorithms  [Buzen  1973,  Reiser  and  Kobayashi  1976,  and  Shum 
1976]  can  be  used  to  calculate  G(N)«  This  makes  the  solution 
technique  efficient. 

In  cases  where  a  subsystem  has  multiple  entries  or 
multiple  exits,  first  the  subsystem  is  converted  to  an 
"equivalent"  single-entry-single-exi t  subsystem.  The 
equivalence  is  defined  as  having  the  same  relative 
throughputs  for  each  centre  in  the  subsystem  in  both  the 
original  and  modified  systems.  This  is  also  characterized 
as  the  transition  matrix  for  the  subsystem  having  a 
"subparallel"  eigenvector  to  the  eigenvector  of  the 
transition  matrix  for  the  whole  system  [Vantilborgh  1978, 
Courtois  1978].  In  single-entry-single-exit  systems, 
"short-circuiting"  is  a  way  of  achieving  this  condition. 


EXAMPLE 


Let  us  consider  an  example  of  multiple-entry-multiple- 
exit  subsystem.  Consider  the  queueing  network  in  Figure 
4.1a.  There  are  three  servers  1,2,  and  3.  Let  Pij,  i, 
j=l, 2, 3,  be  the  probability  of  transition  from  centre  i  to 
centre  j  for  a  job  completing  service  at  centre  i.  Pii  =  0, 
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Figure  4.1b 

Figure  4.1:  Multiple  Entry-Exit  Decomposition 
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i  =  1,2,3.  Suppose  servers  1  and  2  are  combined  together  to 
form  subsystem  A  and  server  3  forms  subsystem  B  (Figure 
4.1b).  Relative  throughputs  for  the  whole  system  are  given 
by : 

el  =  1 

e  2  =  (P32  +  P31*P12)/(P32*P21  +  P32) 
e  3  =  (1  -  P12*P21)/(P32*P21  +  P31) 

For  subsystem  A  the  relative  throughputs  are 
el'  =  1 

e  2  '  =  P12  '/P21 ' 

P12'  and  P21'  are  chosen  such  that  e2  =  e2'  and  0<  P12'  <  1, 
0  <  P21'  <  1.  In  general,  there  exist  many  different 
subsystem  topologies  that  satisfy  these  conditions. 

Proposition  4.1:  If  a  network  satisfies  the  conditions  for 
a  product  form  solution,  (4.1)  is  exact.  Fi(Si)  are  defined 
as  solutions  to  individual  subsytems  (see,  for  example, 
[Baskett  et  al.  1975]). 

Discussion:  If  the  network  satisfies  the  conditions  for 
product  form,  Norton's  theorem  decomposition  yields  an  exact 
solution  [Chandy  et  al.  1975a].  For  the  general  case 
(multi pie -entry -multiple -exit 


situation) , 


the  relative 
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throughputs  for  service  centres  in  a  subsystem  are  the  same 
as  in  the  entire  network.  Thus,  decomposition  gives  an  exact 
result  [Vantilborgh  1978,  Chandy  and  Sauer  1978,  Lazowska 
and  Garner  1979]. 

Solutions  to  individual  subsystems  are  obtained  by 
various  methods  depending  upon  the  arrival  process  of  the 
subsystem.  For  example,  if  each  subsystem  is  composed  of 
one  server  with  a  general  service  time  distribution  and 
arrivals  to  subsystems  are  non-Poisson,  solutions  to 
individual  servers  can  be  obtained  by  solving  them  as  G/G/l 
systems- 

Individual  solutions  to  subsystems,  however,  depend 
upon  the  availability  of  a  solution  technique.  Sometimes 
additional  approximations  are  needed.  For  example,  in  the 
decomposition  discussed  above,  if  no  G/G/l  solver  is 
available,  an  approximation,  such  as  M/G/l,  may  be 
necessary.  At  this  level  we  are  not  confining  ourselves  to 
any  particular  technique.  The  technique  will  depend  upon 
the  situation. 

Another  problem  to  which  we  must  address  ourselves  is 
how  to  estimate  the  arrival  process  at  a  subsystem  (in 
networks  that  do  not  satisfy  the  conditions  for  a  product 
form  solution).  Open  and  closed  networks  are  discussed 
separately.  We  first  discuss  arrival  rates  for  each  centre 
and  then  explain  how  arrival  rates  to  subsystems  can  be 
calculated . 

In  open  systems  absolute  mean  arrival  rates  to 
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individual  service  centres  can  be  calculated  exactly. 
Assume  that  there  are  M  centres,  1 , 2, . . . i , . . , M.  In  steady 
state  the  mean  arrival  rate  is  equal  to  the  mean  departure 
rate  (forced  flow  law  [Denning  and  Buzen  1978]).  Assume 
also  that  Pij  is  the  transition  probability  from  centre  i  to 
j.  Then  external  arrivals  are  assumed  to  be  coming  from 
centre  0  and  the  departures  going  to  centre  M+l. 

Thus,  if  ri  is  the  arrival-  rate  (also  equal  to  the 
departure  rate)  at  the  ith  centre: 


r  i 


L 


Pj i*r j 


(4.2) 


^~ °  for  i  =  l, 2, . . . , M, M+l 


The  rate  of  arrival  from  the  outside,  rO,  and  Pi  j 's  are 
known.  Equation  (4.2)  gives  a  system  of  M  equations  in  M 
unknowns  and  (4.2)  can  be  solved. 

In  closed  systems,  however,  arrival  rates  to  each 
centre  are  known  only  up  to  a  multiplicative  constant  [Buzen 
1971],  For  example,  let  ri  be  the  arrival  rate  at  the  ith 


centre . 


r  i 


M 

Pji*rj 
for  i = 1 , 2 , 


(4.3) 

M 


This  is  not  an  independent  set  of  equations  and  if 
(rl,r 2, . . . ,rM)  satisfies  (4.3),  (K*r 1, K*r 2, . . . , K*rM) 


satisfies  the  above  equations  as  well  for  any  constant  K 
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The  set  of  equations  represented  by  (4.3)  forms  M-l 
independent  equations  in  M  unknowns,  and  has  a  unique 
solution  only  if  some  additional  constraint  is  available. 
The  basic  difference  between  (4.2)  and  (4.3)  is  that  in 

(4.2)  we  know  rO,  and  knowing  rO,  the  ri's  can  be  computed 
exactly.  There  is,  however,  no  such  value  available  in 

(4.3) .  Iterative  techniques  have  been  proposed  to  obtain 
the  absolute  arrival  rates  [Chandy  et  al .  1975b,  Shum  1976, 
Marie  1979],  These  techniques  define  invariants  for  the 
system  and  iterate  the  solution  until  these  invariants  are 
satisfied.  Some  of  the  invariants  commonly  used  are: 

nl  +  n2  +  ....  +  nM  =  N  (4.4) 


£j(((l-Pi(0))/(ri/ci))  -  Kf 
c=l 


0 


(4.5) 


where  ci  is  the  service  rate  at  the  i-th  centre.  The  first 
invariant  denotes  the  fact  that  in  a  closed  system,  the 
total  number  of  customers  is  a  specified  constant.  The 
second  invariant  asserts  that  in  the  steady  state,  the 
forced  flow  law  holds. 

In  the  rest  of  this  section  we  show  that  most  of  the 
approximate  solution  techniques  can  be  viewed  in  the 
framework  of  generalized  product  form.  To  achieve  this  we 
identify  the  subsystems,  techniques  used  to  estimate  the 
arrival  processes,  and  methods  of  solving  the  individual 
subsystems.  The  discussion  here  is  brief  because  these 
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approximations  have  been  described  in  Chapter  2. 


I-  Chandy,  Herzog  and  Woo's  Iteration  [Chandy  et  al -  1975b] 


The  approximation  technique  proposed  by  Chandy  et  al. 
[1975b]  solves  each  service  centre  separately  by  assuming 
the  rest  of  the  system  to  have  a  product  form  solution. 
Thus,  the  arrival  rate  is  dependent  upon  the  number  of 
customers  in  the  subsystem  (the  rest  of  the  system  is 
replaced  by  a  load-dependent  server).  The  arrival  process 
is  a  load-dependent  Poisson  process.  Both  invariants,  (4.4) 
and  (4.5),  are  used  to  achieve  the  desired  level  of 
accuracy,  by  iteration.  At  each  iteration,  the  mean  service 
times  at  servers  are  adjusted  to  reflect  the  effect  of 
higher  moments  of  the  service  time  distribution. 


II.  Shum  and  Buzen's  Extended  Product  Form  [Shum  and  Buzen 
1977] 

The  extended  product  form  technique  is  similar  to 
Chandy  et  al.'s  iteration  technique.  EPF  also  uses  the 
invariants  (4.4)  and  (4.5);  each  service  centre  represents  a 
subsystem.  The  main  difference  is  in  the  way  the  higher 
moments  of  the  service  time  distributions  are  reflected  in 
the  solution.  EPF  does  not  alter  the  mean  service  times  at 
the  centres  to  satisfy  the  invariants  as  discussed  for 
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Chandy  et  al.'s  case.  Instead,  a  search  algorithm  for  the 
absolute  arrival  rate  is  used  to  get  the  best  possible  value 
such  that  error  in  the  invariant  (4.5)  is  minimized.  EPF 
uses  the  convolution  algorithm  for  G(N)  and  therefore 
invariant  (4.4)  is  automatically  satisfied. 

Thus,  each  subsystem  is  solved  as  an  M/G/l/N  system  for 
a  given  arrival  rate.  The  values  of  the  left  hand  side  of 
(4.5)  are  calculated  for  the  entire  range  of  arrival  rates 
and  the  minimum  value  gives  the  desired  solution. 


III.  Marie's  Approximation  Technique  [Marie  1979] 


Marie's  approximation  technique  is  more  general  than 
EPF  and  Chandy  et  al.'s  solution  techniques.  However,  the 
technique  is  very  similar  to  both  of  them.  The  main 
difference  between  EPF  and  this  approximation  is  that  the 
arrival  rates  to  subsystems  are  dependent  upon  the  number  of 
jobs  in  them  in  Marie's  approximation,  but  not  in  EPF. 

According  to  Marie's  description,  a  subsystem  in  this 
approximation  may  have  more  than  one  service  centre.  The 
arrival  process  is  assumed  to  be  Poisson  for  each  subsystem 
with  the  arrival  rate  dependent  upon  the  number  of  jobs  in 
the  subsystem.  To  obtain  the  absolute  arrival  rate,  an 
iterative  procedure  is  used  such  that  invariants  (4.4)  and 
(4.5)  are  satisfied.  The  subsystems  are  solved  using 
numerical  approximation. 
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IV.  Norton's  Theorem  Approximation  [Chandy  et  al.  1975a] 

Norton's  theorem  approximation  (basically 
approximations  suggested  by  Sauer  and  Chandy  [1975]  and 
improved  by  Sevcik  et  al.  [1977])  can  be  thought  of  as  being 
an  extreme  case  of  the  generalized  product  form.  The  whole 
system  forms  the  subsystem;  thus,  the  product  form  (4.1) 
contains  only  one  element. 


V.  Kuhn's  Approximation  [Kuhn  1976] 


Kuhn's  approximation  is  defined  only  for  open  queueing 
networks.  In  open  queueing  networks  the  absolute  arrival 


rate  can  be  computed 

based  upon  the 

rate 

of  arrivals 

from 

outside  (see  equation  (4.2)). 

Kuhn 

characterizes 

the 

transition  processes 

by  the  first  two  moments.  Since 

the 

approximation  is 

defined  only 

for 

open  networks 

no 

normalization  constant  is  needed,  and 

P (nl,n2, . . .nM)  =  Pi (nl) *P2 (n2) *. . . *PM (nM) . 

Each  server  is  a  subsystem,  interactions  are  measured  in 
terms  of  first  two  moments,  and  a  G/G/l  approximation  is 


used  for  each  subsystem. 
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4.3  Generalized  Product  Form:  Two  Realizations 


Having  described  the  generalized  product  form  (GPF)  at 
a  conceptual  level,  we  now  present  two  approximation 
techniques  based  upon  it.  The  approximations  described  here 
are  examples  of  the  generalized  product  form.  In  general, 
the  generalized  product  form  may  be  used  to  develop  a  family 
of  approximations.  A  family  of  these  approximations  can  be 
characterized  by  the  size  of  individual  subsystems,  the 
solution  method  used  to  solve  the  subsystems,  and  the 
methods  used  to  characterize  the  interactions  between  the 
subsystems. 

One  of  the  major  errors  in  using  diffusion 
approximations  for  closed  queueing  networks  is  the  lack  of 
good  estimates  for  the  absolute  arrival  rates  at  the  service 
centres  (see  Chapter  3).  In  our  first  example  of  GPF,  we 
shall  show  that  by  using  the  structure  of  the  generalized 
product  form,  the  diffusion  approximations  can  be 
considerably  improved.  Another  source  of  errors  in  using  a 
product  form  approximation  is  the  lack  of  accurate  solutions 
to  the  individual  subsystems.  In  our  second  example,  we 
shall  use  simulations  to  solve  individual  subsystems  and 
combine  these  solutions  in  the  framework  of  generalized 
product  form. 

This  section  is  divided  into  two  parts:  Part  1 
describes  the  solution  algorithms,  GPF  with  diffusion  and 
GPF  with  simulation,  in  detail,  and  part  2  presents  examples 
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to  investigate  the  accuracy  of  these  new  apr ox imations . 


4-3-1  The  Solution  Algorithm 


We  present  two  solution  algorithms  for  single  class 
queueing  networks  in  detail.  Consider  equation  (4.1).  Si 
represents  the  state  at  subsystem  i.  Diffusion 
approximations,  as  described  in  Chapter  3,  decompose  the 
system  into  subsystems,  each  consisting  of  a  single  service 
centre  in  the  network-  The  state  of  the  ith  subsystem  is 
represented  by  ni,  the  number  of  customers  at  the  ith 
server  - 


If  we  want  to  solve  the  network  exactly,  the 
description  of  a  centre  by  the  number  of  customers  present 
may  not  be  sufficient  to  write  global  balance  equations-  In 
general,  more  than  just  the  number  of  customers  at  the 
centre  is  required  to  get  the  balance  equations-  For 
example,  with  FCFS  and  a  non-exponential  service  time 
distribution,  such  information  as  the  amount  of  service 
attained  by  the  job  currently  being  served  is  needed-  For 
the  purpose  of  calculating  performance  measures  such  as  mean 
queue  lengths  and  server  utilizations,  however,  the 
equilibrium  probability  distribution  for  the  number  of 
customers  at  various  centres  is  Sufficient- 

Substituting  Pi(ni)  for  Fi(Si),  (4-1)  reduces  to: 


P (n 1 , n  2 , « - « , nM) 


( 1/G (N) ) *P1 (nl) *- - -*PM (nM) 


(4-6) 
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Pi  (ni)  is  the  equilibrium  probability  of  ni  customers 
present  at  the  ith  subsystem-  Each  server  is  characterized 
by  a  G/G/l  system.  Two  basic  problems  are: 


1-  Solutions  to  G/G/l  systems 


2.  Estimating  arrival  processes  at  the  subsystems. 

Two  approaches  to  solve  the  first  problem  lead  to  the 
two  examples  of  GPF.  The  first  one  uses  diffusion  and  the 
second  one  uses  simulation.  Several  approaches  to  using 
diffusion  approximations  for  G/G/l  systems  were  discussed  in 
Chapter  2  (for  example,  methods  proposed  by  Kobayashi 
[1974],  Gelenbe  [1975],  and  Yu  [1977]).  In  our  experiments 
with  various  sets  of  examples  and  parameter  values  of  a 
network  we  have  found  Kobayashi's  method  for  solving 
individual  service  centres  to  be  the  most  robust.  (This  has 
also  been  reported  by  Yu  [1977].)  Therefore,  in  GPF  with 
diffusion,  Kobayashi's  [1974]  solutions  are  used  (see 
Chapter  2) . 

In  GPF  with  simulation,  the  general  in ter -ar r ival  time 
distribution  and  the  general  service  time  distribution  are 
approximated  by  the  first  two  moments.  Based  upon  the  value 
of  the  coefficient  of  variation  of  a  distribution, 
exponential,  hyper exponential ,  or  hypoexponen tial 
representations  are  used. 

In  both  approximations,  the  second  problem  is  solved  by 
using  an  iterative  process.  Absolute  arrival  rates  are 
needed  for  each  service  centre  in  the  network.  The  arrival 
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processes  are  characterized  by  the  first  two  moments  of  the 
inter-arrival  time  distribution.  Shum  and  Buzen  [1977] 
estimate  only  the  mean  arrival  rate  at  the  centres.  We  need 
the  second  moments  of  the  inter-arrival  time  distribution  as 
well . 

In  calculating  the  mean  arrival  rate  we  use  the  search 
procedure  described  by  Shum  [1976].  The  entire  search 
interval  is  divided  into  ten  parts.  For  each  boundary  point 
we  calculate  the  left  hand  side  of  the  invariant  (4.5), 
namely , 


(1-Pi  (0)  )/(ri/ci)  -  K)2 


E 


where  K  is  the  value  of  ( (1-Ps (0) )/ (rs/cs) )  for  the  servers 
with  highest  utilization.  Two  intervals  with  a  common 
boundary  at  the  point  where  E  is  minimum  are  further  divided 
into  ten  intervals  and  the  same  procedure  continues.  In  the 
ideal  situation  we  want  E  to  be  0.  In  practice,  however, 
we  stop  the  iteration  once  E  is  less  than  or  equal  to  e, 
some  prespecified  small  value.  For  a  more  detailed 
treatment  of  the  convergence  of  such  an  iteration  see  Shum's 
thesis  [1976]. 

At  each  iteration  for  a  given  arrival  rate  at  the 
service  centre,  the  utilization  is  calculated.  Using  this 
utilization  and  equations  (3.4)  -  (3.6),  the  second  moment 
of  the  arrival  process  is  calculated.  These  values  are  then 
used  to  obtain  diffusion  approximation  solutions  to  each 
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service  centre  in  GPF  with  diffusion  and  for  simulation  in 
GPF  with  simulation.  Equation  (4.6)  is  used  to  obtain 
P  (nl, n 2, . . -nM) .  The  normalization  constant,  G(N), 
guarantees  that  the  total  mean  queue  length  is  N,  thus 
satisfying  the  invariant  (4.4). 


4.3.2  Examples 


We  present  the  set  of  examples  of  queueing  network 
models  discussed  in  Chapter  3,  and  compare  the  accuracy  of 
the  GPF  with  diffusion  and  GPF  with  simulation 
approximations  with  those,  of  EPF  and  modified  diffusion. 
Behaviours  of  these  approximations  are  discussed  in  more 
detail  in  Chapter  5. 

Figures  4.2  and  4.3  present  the  relative  percentage 
errors  in  the  server  utilization  and  the  mean  queue  length 
for  server  1  for  different  numbers  of  service  centres  (M)  in 
the  network.  In  both  these  cases  the  generalized  product 
form  approximations  follow  the  same  error  behaviour  but  with 
smaller  relative  error  than  EPF  and  the  modified  diffusion 
approximation.  Because  of  the  unavailability  of  the  exact 
results  for  M  larger  than  9,  we  are  unable  to  examine  the 
behaviour  of  the  error  curves  for  larger  values.  However, 
to  see  the  impact  of  a  larger  M,  we  use  simulation  results 
for  M  equal  to  10  and  12.  The  errors  compared  to  simulation 
results  are  presented  by  dotted  lines  in  Figures  4.2  and 
4.3.  As  M  increases  all  the  approximations  have  smaller 
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%  Relative  Multiprogramming  Level  (N)  =  3 

Error  in  Mean  Queue  CV  of  the  Service  Time  Distribution 

Length  at  Centre  1  at  the  Central  Server  (CV1)  = 
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Figure  4.3:  Error  in  Mean  Queue  Length  at  CPU  for  Varying  Number  of  Service  Centres  (M) 
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errors  in  utilization  as  well  as  queue  length. 

Figures  4.4  and  4.5  represent  the  behaviour  of  the 
error  for  various  values  of  coefficient  of  variations  of  the 
service  time  distribution  at  the  CPU.  The  error,  again,  is 
reduced  by  using  generalized  product  form  with  diffusion 
instead  of  the  modified  diffusion  approximation.  The 
situation  is  slightly  different  for  GPF  with  diffusion  as 
compared  to  EPF.  For  small  values  of  CV,  GPF  with  diffusion 
reduces  the  error.  However,  for  large  values  of  CV,  EPF 
outperforms  it.  The  EPF  on  the  other  hand,  gives  better 
approximations  for  larger  values  of  CV  of  the  service  time 
d is tr ibution. 

Thus,  although  we  are  characterizing  the  transition 
process  more  accurately  than  EPF,  EPF  is  performing  better 
than  the  proposed  approximation  for  large  values  of  CV  of 
the  service  time  distribution.  The  high  error  for  large 
CV's  is  due  to  poor  accuracy  of  the  diffusion  approximation 
for  G/G/l  in  the  case  of  a  large  CV  [Yu  1977].  This 
hypothesis  is  supported  by  the  error  curves  obtained  for  GPF 
with  simulation.  The  percent  relative  errors,  for  large 
value  of  CV  of  the  service  time  distribution,  using 
generalized  product  form  with  simulation  are  smaller  than 
for  all  other  approximations  for  both  utilization  and  mean 
queue  length.  Using  simulation  results  for  G/G/l  improves 
the  generalized  product  form  considerably  and  performs 
better  than  EPF.  The  technique  using  simulation,  however, 
is  much  more  costly  than  the  one  using  diffusion. 


At  each 


%  Relative 

Error  in  Utilization  Multiprogramming  Level  (N)  = 

at  Centre  1  Number  of  Service  Centres  (M) 


II 


Figure  4.4:  Error  in  CPU  Utilization  for  Varying  CV  at  the  CPU 


%  Relative 

Error  in  Mean  Queue  Multiprogramming  Level  (N)  = 

Length  at  Centre  1  Number  of  Service  Centres  (M) 
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Figure  4.5:  Error  in  Mean  Queue  Length  at  CPU  for  Varying  CV  at  the  CPU 
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step  in  the  search  process  for  absolute  throughput  rate, 
G/G/l  simulations  for  all  the  service  centres  are  required. 
An  alternative  to  simulation  and  diffusion  approximations 
for  each  centre  would  be  numerical  solution  for  a  phase-type 
arrival  and  service  time  distributions  as  suggested  by  Marie 
[1979] . 

The  final  set  of  experiments  is  based  upon  changing  the 
number  of  customers  (multiprogramming  level,  N)  in  the 
system.  The  results  are  presented  in  Figures  4.6  and  4.7. 
Here  again,  using  more  information  than  EPF,  (i.e.,  using 
higher  moments  of  arrival  processes)  improves  the  accuracy. 
The  change  in  N  can  also  be  interpreted  as  change  in  load. 

4.4  Summary  and  Comments 

In  this  chapter,  we  have  presented  a  general 
approximate  solution  technique  for  queueing  networks.  It 
was  shown  that  most  of  the  existing  approximations  can  be 
viewed  in  the  framework  of  the  general  approximation,  GPF. 
As  an  example  of  the  usefulness  of  the  concept  of  defining  a 
general  technique,  two  new  approximations  were  presented. 
These  techniques  are  very  closely  related  to  EPF  and 
diffusion  approximation  techniques.  With  the  help  of  a  set 
of  examples,  it  was  shown  that  the  new  approximation  has 
better  accuracy  than  both  EPF  and  diffusion  in  most  cases 
over  a  range  of  parameter  values. 


°l  Relative  Number  of  Service  Centres  (M)  =  4 

Error  in  Utilization  C V  of  the  Service  Time  Distribution 

at  Centre  1  at  the  Central  Server  (CV1)  =  4 
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Figure  4.6:  Error  in  CPU  Utilization  for  Varying  Multiprogramming  Level  (N) 


%  Relative  Number  of  Service  Centres  (M)  =  4 

Error  in  Mean  Queue  C V  of  the  Service  Time  Distribution 

Length  at  Centre  1  at  the  Central  Server  (CV1)  =  4 
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Figure  4.7  :  Error  in  Mean  Queue  Length  at  CPU  for  Varying  Multiprogramming  Level  (N) 
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Two  questions,  namely,  how  much  extra  computational 
cost  is  incurred  to  get  the  extra  accuracy  for  GPF  over 
other  approximations  and  why  is  GPF  with  simulation  any 
better  than  simulating  the  entire  network,  need  further 
discussion-  GPF  with  diffusion  requires  no  more  computation 
than  that  required  by  EPF-  The  only  difference  in  cost 
between  modified  diffusion  and  GPF  with  diffusion  is  the 
number  of  iterations  required  for  GPF  with  diffusion.  In 
examples  presented  in  this  chapter,  the. average  number  of 
iterations  was  5-  In  terms  of  absolute  computational  cost 
five  iterations  do  not  make  any  significant  difference-  For 
example,  both  approximations  were  executed  using  the  high 
speed  job  stream  (HSJS)  of  the  University  of  Toronto 
Computer  Centre-  (Note  that  the  maximum  CPU  execution  time 
allowed  in  HSJS  was  3  seconds  on  an  IBM  370/168- )  The 
advantage  of  GPF  with  simulation  over  simulating  the  entire 
model  is  reduced  cost.  The  time  required  to  attain  a  given 
level  of  accuracy  by  simulating  the  entire  network  is  many 


times  more  than  that  required  by  simulating 
centres  times  the  number  of  centres- 


the  individual 
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Chapter  5 

ERROR  BEHAVIOUR:  AN  EXPERIMENTAL  STUDY 

This  chapter  presents  a  comparative  study  of  various 
approximate  solution  techniques-  The  performance  of  the 
approximations  is  compared  based  upon  three  basic 
parameters : 

1-  coefficient  of  variation  of  service  time  distribution 
(CV) 

2.  number  of  customers (N) 

3.  number  of  servers  (M) 

This  list  is  not  exhaustive.  However,  these  parameters 
reflect  the  basic  changes  one  would  like  to  make  in  a 
network-  For  example,  changing  M  and  N  have  compounding 
effects-  If  M  is  increased  for  a  fixed  N,  the  load  at 
various  servers  will  generally  be  decreased.  Similarly,  by 
holding  M  fixed  and  decreasing  N,  a  similar  effect  can  be 
achieved.  Also,  by  varying  M  and  N,  a  third  effect,  namely 
loading,  can  be  reflected,  as  well- 

One  of  the  other  properties  of  the  network,  namely 
balanced  vs-  non-balanced ,  can  also  be  reflected  by  the 
parameters  mentioned  above-  By  increasing  the  number  of 
centres,  M,  the  CPU  saturates  much  faster  than  other  servers 


-112- 


in  the  system-  Thus,  an  approximately  balanced  system  for  a 
smaller  M  yields  an  unbalanced  system  for  a  large  M« 

The  set  of  networks  to  be  considered  should  include 
both  realistic  and  extreme  cases  -  realistic  in  the  sense 
that  case  studies  of  the  queueing  network  models  of  actual 
computer  systems  resemble  them,  extreme  in  that  all  boundary 
conditions  are  included  to  demonstrate  robustness. 

In  most  of  our  experiments  the  topology  of  the  queueing 
networks  has  been  restricted  to  a  closed  central  server 
model-  Central  server  models  are  widely  used  in  case 
studies  and  they  reflect  the  physical  behaviour  of  many 
systems  - 

As  discussed  in  Chapter  4,  closed  networks  are  more 
difficult  to  solve  than  open  networks.  We  have  ignored  open 
networks  in  our  comparisons-  Also,  most  of  the  approximate 
solution  techniques  compared  are  only  defined  for  closed 
networks,  for  example,  EPF,  Norton's  theorem  decomposition 
(Sauer  and  Chandy  [1975]  and  Sevcik  et  al-  [1977])-  The 
following  approximations  are  compared: 

1-  Norton's  theorem  by  SLTZ  [Sevcik  et  al .  1977] 

2-  Modified  diffusion  approximations  [Chapter  3] 

3.  EPF  [Shum  1976] 

4-  GPF  with  diffusion  [Chapter  4] 


5. 


GPF  with  simulation  [Chapter  4] 
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In  our  experimentation  we  use  the  same  central  server 
model  as  presented  in  Chapter  3-  The  basic  parameters  are, 
however,  varied  in  these  experiments.  Changes  made  are 
discussed  where  applicable. 

This  chapter  is  divided  into  four  sections.  The  first 
three  sections  present  the  effect  of  the  three  parameters 
mentioned  above  (CV  of  the  service  time  distribution,  N,  and 
M)  on  the  error  in  estimating  the  utilizations  and  mean 
queue  lengths.  The  presentation  of  results  in  these 
sections  are  in  similar  format:  the  error  behaviours 
followed  by  explanations  of  the  observed  behaviours. 
Finally,  Section  4  presents  a  summary  of  the  results 
obtained. 


5.1  Service  Time  Distribution 

Service  time  distributions  have  received  more  attention 
in  evaluating  approximate  solution  techniques  than  many 
other  parameters  [Sauer  1975,  Shum  1976,  Gelenbe  1975, 
Reiser  and  Kobayashi  1975].  For  networks  to  have  exact 
product  form  solutions,  either  exponential  service  time 
distributions  at  service  centres  are  assumed,  or  station 
balance  [Chandy  et  al.  1977]  service  disciplines  (such  as 
PS,  PRLCFS,  NQ) ,  which  eliminate  all  the  effects  of 
distributional  form  other  than  the  mean,  are  assumed.  For 
FCFS ,  an  exponential  service  time  distribution  is  needed. 
This,  however,  is  contrary  to  many  realistic  situations  for 
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computer  systems  [Sauer  1975,  Baskett  1971,  Lester  19731- 

In  this  section  we  present  the  error  behaviours  of 
various  approximations  for  a  wide  range  of  CV's  of  service 
time  distributions.  To  obtain  an  exact  steady-state 
solution  for  the  network,  a  Cox-type  (phase-type) 
representation  for  the  servers  is  assumed.  This  section  is 
divided  into  three  parts.  Part  1  presents  the  error  curves, 
part  2  explains  the  results,  and,  finally,  part  3  comments 
on  simulating  servers  with  high  CV's. 


5.1.1  Experimental  Results 


Figure  5.1  presents  the  error  behaviours  of  various 
approximations  in  estimating  the  utilization  of  the  central 
server  for  a  range  of  CV  values.  The  range  presented  here 
is  based  upon  two  considerations:  the  lower  bound  (.25)  was 
due  to  limitations  on  the  state-space  size  to  obtain  exact 
solutions  by  QSOLVE;  the  upper  bound  (225)  basically  shows 
the  trend  and  no  further  insights  were  obtained  by 
increasing  this  range.  Figure  5.2  presents  similar  results 
for  mean  queue  length  at  the  central  server  over  the  same 
range  of  CV. 

The  general  behaviour  of  the  two  curves  is  similar  for 
all  the  approximations:  the  actual  errors  in  mean  queue 
lengths  are  substantially  higher  than  those  for  utilizations 
in  all  approximations.  The  behaviour  of  various 


approximations,  however, 


is  not  the  same  for  the  entire 
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Figure  5.1:  Error  in  CPU  Utilization  for  Varying  CV  at  the  CPU 
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Figure  5.2:  Error  in  Mean  Queue  Length  at  CPU  for  Varying  CV  at  the  CPU 
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range  - 

For  CV>1 ,  as  we  increase  the  CV,  the  errors  tend  to 
increase  in  all  cases.  However,  for  EPF  and  Norton's 
theorem  by  SLTZ  the  error  levels  off  and  then  decreases  for 
CV  beyond  100.  Similar  observations  for  high  CV's  have  been 
reported  by  Tripathi  and  Levy  [1976]  for  Norton's  theorem 
(Sauer  and  Chandy's  case)  and  by  Shum  [1976]  for  EPF.  The 
leveling-off  point  varies  from  network  to  network  based  upon 
other  parameters,  but  the  same  trend  is  obtained  for  all 
networks.  For  the  approximations  based  upon  diffusion, 
namely  GPF  with  diffusion  and  modified  diffusion 
approximations,  the  error  increases  for  higher  values  of  CV 
and  no  decreasing  trend  was  found  in  any  of  the  experiments. 
The  error  for  both  utilization  and  mean  queue  length  using 
GPF  with  simulation  is  considerably  smaller  than  those 
arising  from  all  other  approximations. 

For  CV<1 ,  as  the  CV  decreases,  the  errors  tend  to 
increase  slightly  although  the  magnitude  is  much  smaller 
than  that  for  CV>1  and  differs  in  various  approximations. 
No  anomalies  were  observed  in  this  case. 

Overall,  GPF  with  simulation  has  the  smallest  error  and 
GPF  with  diffusion  performs  well  for  moderate  values  of  the 
CV  of  the  service  time  distribution  (less  than  100).  For 
larger  values  of  CV,  however,  GPF  with  diffusion  has  large 
errors  in  both  utilization  and  mean  queue  length.  The 
behaviour  observed  in  these  experiments  is  explained  in  the 


next  section. 
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5-1-2  Explaining  the  Error  Curves 


Increasing  CV  above  one  or  decreasing  it  below  one  for 
the  FCFS  service  discipline  in  effect  deviates  from  the 
conditions  the  network  must  satisfy  to  yield  a  product  form 
solution.  Intuitively,  the  error  curves  should  reflect  this 
property.  That  is,  the  errors  in  performance  measures  may 
be  increased  in  those  cases.  As  it  turns  out,  for  all  the 
approximations,  going  away  from  CV  equal  to  one,  the  errors 
are  increased  except  in  EPF  and  in  Norton's  theorem  by  SLTZ. 
The  behaviour  of  GPF  with  simulation  does  not  show  any 
significant  change  due  to  variation  in  CV  of  the  service 
time  distribution. 

The  behaviour  of  GPF  with  diffusion  needs  more 
explanation.  GPF  with  diffusion  outperforms  all  other 
non-GPF  approximations  for  moderate  values  of  CV  but 
performs  worse  than  EPF  or  Norton's  theorem  by  SLTZ  for 
large  CV's.  One  of  the  basic  reasons  for  this  is  the  poor 
performance  of  the  diffusion  approximation  of  G/G/l  with  a 
high  CV  of  the  service  time  distribution.  (For  a  more 
detailed  explanation  see  Chapter  4.)  The  other  reason  will 
be  clear  once  we  discuss  the  behaviours  of  EPF  and  Norton's 
theorem  by  SLTZ  for  large  CV's. 

To  explain  the  error  behaviours  in  EPF  and  Norton's 
theorem  by  SLTZ  approximations  for  large  CV's,  let  us  first 
recall  the  representations  for  distributions  with  a  large  CV 
(>1).  As  discussed  in  Chapter  2,  if  M  is  the  mean  and  CV  is 
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the  coefficient  of  variation  of  a  distribution  with  a  CV 
greater  than  1,  then  the  representation  in  Figure  5.3  is 
used . 


Figure  5.3 

P,  Ml,  and  M2  can  be  determined  by  solving  the  following  set 
of  equations. 


P*M1  +  ( 1— P) *M  2  =  M 


(5.1) 

2* [P*M1*M1  +  (1-P) *M2*M2]/(M*M)  =  CV+1 


There  are  an  infinite  number  of  solutions  to  the  system 
of  equations  given  by  (5.1),  because  there  are  two  equations 
in  three  unknowns.  It  has  been  shown  that  the  selection  of 
a  particular  solution  to  (5.1)  can  affect  the  values  of 
performance  measures  [Price  1976,  Lazowska  1977].  For 
example,  Lazowska  [1977]  has  shown  that  two  different 
solutions  to  (5.1)  can  produce  results  for  server 
utilizations  that  differ  by  more  than  20%. 

To  solve  (5.1),  an  extra  condition  must  be  introduced. 
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For  example,  Sauer  [1975]  suggested  that  the  load  at  both 
the  servers  should  be  the  same  -  That  is, 

P*M 1  =  (1-P) *M2  (5.2) 

Shum  [1976],  based  upon  Price's  [1976]  work,  recommended 
minimizing  the  difference  between  the  means  of  two 
exponential  distributions.  That  is, 

Minimize [abs (Ml  -  M2)]  (5.3) 

Both  (5.2)  and  (5.3)  lead  to  the  same  solution 

P  =  (1/2) * [1  -  sqrt( (CV-1) / (CV+1) ) ] 

Ml  =  M/(2*P)  (5.4) 

M2  =  M/  (2* (1-P) ) 

For  a  fixed  mean  and  a  varying  CV  of  a  hyper exponential 
distribution,  the  resulting  values  of  P  are  given  in  Table 
5.1.  Approximate  values  for  maximum  P  (Pmax)  are  also  given 
in  the  table.  Note  that  as  CV  increases,  P  very  quickly 
approaches  zero.  Furthermore,  the  range  in  which  P  can  vary 
decreases  with  increasing  CV  [Lazowska  1977]. 


CV 

P 

Pmax 

4 

.  1127 

.  4000 

9 

.  0527 

.  2836 

16 

.  0303 

.  1354 

64 

.  0077 

.  0834 

100 

.  0049 

.  0095 

400 

.  0012 

.  0052 

900 

.  0005 

.  0017 

2,  500 

.  0002 

.  0005 

Table  5.1: 

P  for  high  CV 

Thus,  with  the 

two  stage  representation,  for  a  very 

large 

CV,  the  choice  of  P  is  limited  and  (5.4)  gives  values 

close 

to  zero.  In 

other  words,  for 

a  centre  with  a  service 

time 

distr ibution 

having  a  large 

CV,  the  vast  majority  of 

jobs 

experience 

the  same  exponential  service  time 

distribution  (namely  the  one  chosen  with  1-P  probability). 
Only  for  very  few  jobs  is  the  service  time  drawn  from  the 
distribution  with  a  large  mean.  But,  since  there  are  only  N 
customers  only  a  "limited  damage"  [Buzen  1977,  Lazowska 
1977]  to  other  jobs  can  be  done.  The  damage  due  a  long  job 
can  occur  only  infrequently  and  at  most  N-l  jobs  can  be 
affected  at  any  of  those  instances.  Thus,  the  server  with  a 
service  time  distribution  having  a  very  high  CV  very  closely 
resembles  a  server  with  an  exponential  service  time 


dis tr ibution. 


This  suggests  that  any  technique  that  uses 
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such  a  representation  for  non-exponential  servers  and 
performs  well  (less  error)  for  exponential  cases  will 
perform  well  for  very  high  CV's. 

Finally,  note  that  no  such  representations  are  used  in 
diffusion  approximations-  The  following  conjecture  can  be 
made:  "If  an  approximation  technique  makes  use  of  the 
representational  form  of  the  hyperexponential  distribution 
mentioned  above  and  has  small  error  for  the  exponential 
distribution  the  error  will  be  small  for  very  large  CV's." 
The  results  presented  here  and  observations  made  by  Shum 
[1976]  support  this  conjecture. 

5-1-3  Comments  on  Simulating  Servers  with  High  CV 

In  many  simulation  studies  of  queueing  systems,  servers 
with  specified  first  two  moments  of  the  service  time 
distributions  are  represented  by  a  Cox-type  server-  For 
example,  if  the  CV  of  the  service  time  distribution  is 
greater  than  one  a  two-stage  hyperexponen t ial  representation 
is  used  (Figure  5.3). 

Based  upon  the  results  presented  in  the  last  section, 
we  show  here  that,  in  simulating  a  server  with  high  CV,  a 
very  large  number  of  transactions  are  needed-  Let  us 
consider  the  server  shown  in  Figure  5-3-  Assume  that  we 
want  to  simulate  a  server  such  that  the  observed  value  of  P 
is  within  10%  of  what  we  want  it  to  be  theoretically-  (The 
importance  of  such  a  restriction  is  discussed  towards  the 
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end  of  this  section.)  For  example,  if  P  *  .1,  the  observed  • 
value  of  P  should  lie  within  (.1  +.01)  with  a  very  high 

probability.  We  want  to  estimate  the  number  of  transactions 
needed  to  satisfy  this. 

The  width  of  the  confidence  interval  at  a  95%  level  of 
confidence  is  [Fishman  1973]: 

W  *  2* (1.96) *sqr t (P* (1-P) /n) 

Thus  an  estimate  of  n,  n',  is 

n '  =  (3.92) * (3.92) *P* (1-P)/(W*W) 

=  (15-3664) *P*(1-P)/(W*W) 

Since  the  observed  value  of  P  should  be  within  10%  of  its 
actual  value,  the  confidence  interval  has  width  equal  to  20% 
of  the  value  of  P,  i.e.,  W  =  P/5.  Therefore, 

n'  =  (384. 16) * (1-P)/P 

For  a  given  value  of  the  CV,  one  can  calculate  P  (see  Table 
5.1)  and  given  a  P,  n '  can  be  obtained  by  the  above  formula. 
The  values  of  n'  for  various  values  of  CV's  are  given  in 
Table  5-2.  Note  that  decreasing  P  below  .1  increases  the 
required  number  of  observations.  Halving  P  more  than 


doubles  n'. 
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cv 

n  ' 

4 

3,  000 

9 

7,  000 

16 

13, 000 

64 

50, 000 

100 

77, 000 

400 

308, 000 

900 

700, 000 

2,  500 

2, 000, 000 

Table  5.2:  n'  for  high  CV's 


These  values  of  n'  are  lower  bounds  because  we  have 
assumed  all  the  regularity  conditions  needed  for  the  central 
limit  theorem  [Fishman  1973].  In  general,  more  transactions 
will  be  needed.  It  should  also  be  noted  that  if  we  want  to 
reduce  the  width  of  the  confidence  interval,  the  number  of 
observations  needed  are  proportional  to  the  square  of  the 
reduction  factor.  For  example,  to  have  a  confidence 
interval  with  a  width  equal  to  half  of  the  above  confidence 
interval,  4*n'  observations  have  to  be  taken.  In  a  queueing 
network  all  the  transactions  do  not  go  through  all  the 
servers.  Only  a  fraction,  depending  upon  routing  frequency, 
pass  through  a  given  server.  Thus,  an  even  larger  number  of 
transactions  is  needed. 

The  error  in  estimating  P  can  directly  affect  the 
accuracy  of  performance  measures  (e.g«,  utilization,  mean 
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queue  length).  We  present  a  simple  M/G/l  simulation  example 
to  demonstrate  this  point.  Let  the  service  time 
distribution  have  mean  equal  to  10ms  and  CV  equal  to  900 
{(standard  deviation/mean) =30} .  Assume  that  the  mean 
inter-arrival  time  is  15ms.  Therefore,  the  utilization,  U, 
is  equal  to  2/3.  The  mean  queue  length,  E  [Q] ,  can  be 
calculated  using  the  Pollaczek-Khinchin  equation  [Kleinrock 
1976] : 

E [Q]  =  U  +  (U*U) *{ (l+CV)/(2* (1-U) ) } 

=  601.33 

Now  let  us  consider  a  simulation  experiment  for  this 
M/G/l  queue.  It  is  common  to  represent  a  general  service 
time  distribution  with  CV  greater  than  one  by  a  two  stage 
hyperexponential  distribution  [Lavenberg  et  al.  1979]. 
Using  equation  (5.4)  the  parameters  of  this  hyperexponential 
distribution  can  be  obtained: 


0005, 

Ml  = 

10000,  and  M2 

=  5.0025 

is 

very 

small,  unless 

a  large  number  of 

transactions  is  simulated  a  high  error  in  P  is  not  unlikely. 
Suppose  that  the  observed  value  for  P  in  a  simulation  is 
.0004  (a  twenty  percent  error).  The  mean  and  the  CV  of  the 
simulated  server  can  be  calculated: 

Mean  =  9.0004 


CV  =  987.  18 
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Therefore,  the  estimated  value  for  the  utilization  is  .6000. 
The  mean  queue  length  is  given  by: 

E [Q]  =  445.28 

Thus,  the  error  in  estimating  the  mean  queue  length  is  about 
26%.  A  twenty  percent  error  in  P  has  resulted  in  a  more  than 
twenty  five  percent  error  in  the  mean  queue  length.  Note 
that  a  twenty  percent  error  is  not  unlikely  if  less  than 
700,000  transactions  are  simulated. 

5.2  Number  of  Customers 

0 

In  our  experiments  we  vary  the  number  of  customers 
(also  referred  to  as  multiprogramming  level),  N,  from  1  to 
8.  The  upper  limit  on  N  is  due  to  the  limitations  of  the 
solution  system,  QSOLVE  [Levy  1977],  we  use.  To  avoid  the 
intermingling  of  the  errors  due  to  approximation  and  to 
simulation,  we  confined  ourselves  to  QSOLVE.  The  general 
trend,  however,  is  clear  from  the  network  parameters 
considered. 


5.2.1  Experimental  Results 


Figure  5.4  represents  the  error  behaviour 
utilization  in  five  approximations  -  EPF,  modified 
approximation,  Norton's  theorem  approximation  by 


for  CPU 
diffusion 
SLTZ,  GPF 


%  Relative  Number  of  Service  Centres  (M)  =  4 

Error  in  Utilization  CV  of  the  Service  Time  Distribution 

at  Centre  1  at  the  Central  Server  (CV1)  =  4 
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Figure  5.4:  Error  in  CPU  Utilization  for  Varying  Multiprogramming  Level  (N) 
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with  diffusion,  and  GPF  with  simulation-  All  these 
approximations  are  applied  for  the  same  queueing  network 
with  varying  levels  of  multiprogramming  (N). 

For  all  approximations  the  error  in  the  utilization  of 
server  1  is  small  (within  10%)  for  low  and  high  values  of  N« 
The  critical  values  of  N,  where  the  errors  are  large,  lie  in 
the  middle.  In  the  example  presented,  these  values  are 
between  3-6.  However,  as  discussed  in  the  following,  this 
interval  varies  for  different  cases. 

It  is  easy  to  see  that  if  there  is  only  one  customer  in 
the  system  no  approximations  are  needed.  We  can  assume  any 
scheduling  discipline  for  simplification.  For  example,  PS 
and  FCFS  are  the  same  for  single  customer  models  and 
therefore  we  can  solve  the  model  exactly  for  any  service 
time  distribution.  In  other  words,  there  is  no  queueing  in 
this  situation. 

If  there  are  many  more  customers  than  the  number  of 
service  centres  in  the  network,  the  bottleneck  service 
centre  in  the  network  may  be  close  to  saturation. 
Therefore,  the  state  of  the  network  is  very  highly  governed 
by  this  one  server. 

The  corresponding  error  curves  for  mean  queue  length  at 
the  central  server  are  given  in  Figure  5.5.  Modified 
diffusion  approximation  has  the  largest  error  in  the  mean 


queue  length.  In  general,  the  mean  queue  length  shows 
larger  error  in  all  cases  except  GPF  with  simulation  where 
both  errors  are  small.  The  general  trend  is  the  same  except 
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Figure  5.5:  Error  in  Mean  Queue  Length  at  CPU  for  Varying  Multiprogramming  Level  (N) 
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larger  errors  than  those  for  utilization  are  encountered. 


5.2.2  Asymptotic  Behaviour 

The  unimodal  behaviours  of  errors  in  various 
approximations  for  varying  numbers  of  customers  can  be 
explained  using  asymptotic  arguments.  Wong  [1975]  gives  a 
derivation  for  the  value  of  number  of  customers,  Nmax, 
beyond  which  queueing  is  inevitable.  Nmax  is  also  known  as 
the  system  saturation  point  [Denning  and  Buzen  1978]. 

Nmax  =  (U1+U2+. . . +UM) /Us 

where  Ui  is  the  relative  utilization,  Xi/Mi  (Xi  is  a 
solution  to  (4.2)  and  Mi  is  the  mean  service  time  at  server 
i),  for  the  ith  server  and  Us  is  the  relative  utilization 
for  the  most  highly  utilized  server  in  the  network. 

This  result  is  derived  using  the  asymptotic  bound 
argument.  If  N  is  equal  to  one,  we  know  that  the  exact 
values  of  the  server  utilizations  are  in  the  same  proportion 
as  relative  utilizations  and  sum  to  one.  As  N  increases  the 
upper  bound  increases  linearly.  For  a  saturated  server  s, 
Us  increases  the  fastest.  However,  it  cannot  be  more  than 
one.  If  Us(l)  is  the  utilization  of  the  saturated  server  s 
when  there  is  one  customer  in  the  network,  then 
Nmax=  ( 1/Us ( 1 ) ) .  This  yields  the  above  equation. 

For  N  smaller  than  Nmax,  there  is  less  queueing  at  the 
service  centres  and  therefore,  there  may  be  less  error  in 
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estimating  the  performance  measures.  For  N  larger  than 
Nmax,  the  queueing  is  increased  and  the  bottleneck  server  is 
likely  to  be  saturated.  Therefore,  there  may  be  less  error 
in  estimating  the  performance  measures  for  larger  values  of 
N,  as  well.  For  our  experimental  parameters  (Ul=l,  U2=.8, 
U3=«  9,  and  U4  =  l) : 

Nmax  =3.7 

In  Figure  5.4  all  the  approximations  have  maximum  error  for 
an  N  larger  than  3.7.  This  supports  the  view  that  for  very 
little  queueing  (N<<Nmax)  the  error  in  any  reasonable 
approximation  should  be  small  but  for  N  close  to  Nmax 
approximations  may  have  their  highest  errors. 

5.3  Number  of  Service  Centers 

Keeping  the  topological  structure  of  the  queueing 
network  as  a  central  server  model,  the  number  of  servers  can 
be  increased  by  increasing  the  number  of  I/O  centres.  Also, 
to  avoid  the  effect  due  to  saturating  the  CPU,  all  the 
branching  probabilities  are  readjusted;  all  the  I/O  centres 
experience  lighter  loads.  For  maintaining  approximately  the 
same  load  at  the  I/O  servers,  N  has  to  be  increased.  This, 
however,  leads  to  a  saturated  CPU.  In  our  experiments  we 
have  confined  ourselves  to  only  changing  M«  All  the  new 
added  servers  have  branching  probabilities  equal  to  P2,  the 
branching  probability  to  the  second  server,  and  the  same 
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service  time  distribution  as  server  2.  The  ratio  between  P2 
and  P4  is  maintained  (see  Figure  3-9,  page  81  ). 

This  section  is  divided  into  two  parts:  The  first  part 
presents  the  results  obtained  by  the  experimental  studies, 
and  the  second  part  gives  explanations,  theoretical  and 
intuitive,  for  those  results. 


5-3-1  Experimental  Results 


Figures  5.6  and  5-7  present  the  error  curves  for 
utilization  and  mean  queue  lengths  at  the  central  server  for 
various  approximations-  Since  we  are  considering  queueing 
networks,  there  are  at  least  two  service  centres-  The 
results  are  reported  for  two  to  nine  service  centres.  For 
larger  M  we  found  no  new  trend  and  no  more  insight- 

The  error  behaviour,  in  both  figures,  is  unimodal-  GPF 
with  diffusion  performs  better  than  non-GPF  for  all  values 
of  M  except  for  M  equal  to  two  where  Norton's  theorem 
approximation  by  SLTZ  has  no  error.  This  is  because  a  two 
server  system  is  solved  exactly  [Sauer  1975].  Typical 
errors  in  utilizations  are  between  1%  to  10%,  whereas  errors 
in  mean  queue  lengths  are  of  the  order  5%  to  20%. 


5-3-2  Explaining  Error  Curves 


Before  we  explain  the  behaviour  of  errors  we  prove  two 
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Figure  5.6:  Error  in  CPU  Utilization  for  Varying  Number  of  Service  Centres  (M) 
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Figure  5.7:  Error  in  Mean  Queue  Length  at  CPU  for  Varying  Number  of  Service  Centres  (M) 
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lemmas-  Given  n  renewal  processes,  let  a  new  process  be. 
formed  by  combining  all  their  renewal  epochs  into  one 
sequence-  We  assume  that  the  new  process  is  a  renewal 
process  although  that  is  generally  not  so.  Let  W  be  the 
waiting  time  for  the  first  renewal  epoch.  We  show  that  W  is 
approximately  exponentially  distributed,  and  therefore,  the 
combined  process  is  close  to  a  Poisson  process-  Let  Mi  and 
Fi  be  the  mean  and  the  distribution  function  for  the 
inter-arrival  time  of  the  ith  renewal  process  and  M  and  F  be 
the  mean  and  the  distribution  function  for  the  inter-event 
time  of  the  combined  process. 

Lemma  5.1:  Under  the  assumption  that  each  renewal  process  is 
rare  (so  that  the  cumulative  effect,  the  resulting  process, 
is  due  to  many  small  causes) ,  as  n  increases  F  approaches 
exponential- 

Proof:  See  Appendix  C. 

Let  us  now  consider  the  reverse  case,  i.e«,  a  renewal 
process  is  split  into  n  processes-  Let  Pi  be  the 
probability  (Bernoulli  trials)  of  transition  to  the  ith 
branch,  and  M  and  CV  be  the  mean  and  coefficient  of 
variation  of  the  inter-event  time  of  the  renewal 
process.  Let  Mi  and  CVi  represent  corresponding  mean  and  CV 
of  the  ith  resulting  process. 
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Lemma  5.2:  If  for  every  i  as  n  increases  Pi  decreases 
uniformly  then  CVi  approaches  1  as  n  approaches  infinity. 

Proof:  See  Appendix  C. 


Note  that  we  have  been  characterizing  distributions 
only  by  their  first  two  moments  and  that  if  CV  =  1,  we  have 
assumed  exponential.  Thus,  in  our  representation,  Lemma  5.2 
implies  that  the  resulting  processes  approach  Poisson-  In 
general,  however,  CV=1  does  not  necessarily  imply  the 
exponential  distribution. 

For  large  values  of  M  in  the  network  considered,  the 
decreased  error  in  performance  measures  can  be  explained  by 
the  above  two  lemmas.  EPF  assumes  that  the  arrival  process 
at  each  centre  is  Poisson.  As  M  increases  the  CV  of  the 
inter-arrival  time  distribution  at  each  I/O  server 
approaches  1.  (We  are  using  Lemma  5.2  and  the  fact  that  we 
are  considering  a  central  server  model.)  The  CV  of  the 
inter-arrival  time  distribution  at  the  CPU,  the  central 
server,  also  approaches  1  (by  Lemma  5.1).  Thus,  for  all  the 
servers  in  the  network  the  arrival  processes  are  "nearly 
Poisson".  This  satisfies  the  assumption  made  in  EPF  and 
results  in  better  approximations. 

To  explain  the  behaviour  of  Norton's  theorem 
approximation  by  SLTZ,  let  us  consider  the  I/O  subsystem  in 
isolation.  For  any  arrival  process  at  the  I/O  subsystem,  if 
the  number  of  I/O  servers  is  large  (i.e«,  M  is  large),  the 
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CV  of  the  inter-departure  time  distribution  from  the 
subsystem  is  also  approximately  1  (by  Lemma  5.1).  Thus,  the 
subsystem  satisfies  the  M  =>  M  property  and  a  smaller  error 
in  performance  measures  by  using  the  approximation  is 
expected  [Sevcik  et  al.  1977]. 

Also,  increasing  M  implies  a  higher  load  at  the  central 
server  and  thus  smaller  errors  for  diffusion  approximations 
(see  [Yu  1977]).  The  errors  are  smaller  for  I/O  devices  as 
well.  One  of  the  main  reasons  is  that  all  I/O  devices  and 
the  CPU  experience  arrival  processes  with  CV  approximately  1 
for  large  M  (Lemmas  5.1  and  5.2)  and  diffusion 
approximations  work  better  for  M/G/l  than  G/G/l. 

5.4  Summary 

We  have  presented  results  of  experiments  performed  to 
test  the  accuracy  of  various  approximations  for  a  range  of 
parametric  values.  Intuitive  and  theoretical  (where 
possible)  explanations  of  the  observed  behaviours  have 
been  presented. 

The  experimental  results  presented  in  this  chapter 
support  the  claim  made  in  the  unifying  hypothesis  presented 
in  Chapter  3.  Errors  in  approximations  were  classified  in 
two  categories  and  it  was  proposed  that  an  approximation 
should  reduce  both  of  these  errors  (local  and  global)  to 
have  a  small  overall  error  in  the  solution.  Table  5.3  gives 

local  solutions  to  the  individual 


a  summary  of  how 
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subsystems  are  obtained  and  at  what  level 
among  the  subsystems  are  measured  for 
approximations  compared. 


the 

all 


interactions 
the  five 


Table  5.3 


Approximations 

Local  solutions 

Interactions 

GPF  with 

Simulation 

Simulation 

Mean  and  CV  of  the 
service  time  dist. 

Mean  and  CV 
by  iterative 
process 

GPF  with 

Diffusion 

Diffusion 

Mean  and  CV  of  the 
service  time  dist. 

Mean  and  CV 
by  iterative 
process 

EPF 

Analytic 

Mean  and  CV  of  the 
service  time  dist. 

Mean  by 

iterative 

process 

Modified 

Diffusion 

Diffusion 

Mean  and  CV  of  the 
service  time  dist. 

Mean  by 
assuming 
product  form 
and  CV  by 
equation 

Norton ' s 

theorem  by 

Global  balance 

Mean  by 

SLTZ 

Mean  and  CV  of  the 
service  time  dist. 

assuming 
product  form 
and  CV  by 
iterative 
process 

Based  on  the 
rather  hard  to  gi 
the  various  approx 
are  highly  corre 


information  in  the  above 
ve  any  partial  ordering  on 
imations  because  local  and 
lated.  However,  Table  5. 


table,  it  is 
the  accuracy  of 
global  errors 
3  can  be  used  to 


explain  the  experimental  results  reported  in  this  chapter 
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Chapter  6 


SUMMARY  AND  DIRECTIONS  FOR  FUTURE  RESEARCH 


6-1  Summary 


The  main  results  of  this  thesis  have  been  summarized  in 
the  final  sections  of  Chapters  3,  4,  and  5.  We  shall 

briefly  reiterate  them  here¬ 
in  Chapter  3,  a  unified  view  of  various  approximate 
solution  techniques  is  presented,  based  upon  the  possible 
sources  of  errors  in  approximations.  Most  of  the 
approximation  techniques  are  characterized  as  decomposition 
methods  and  two  kinds  of  errors,  local  and  global,  are 
discussed-  This  characterization  is  applied  to  diffusion 
approximations-  Errors  in  diffusion  approximations  can  be 
reduced  by  reducing  the  global  error.  An  improved 

characterization  of  transition  processes  is  presented. 
Examples  of  both  open  and  closed  queueing  network  models 
show  the  improvement  in  the  existing  diffusion 

approximations  by  using  the  derived  results  for  transition 
process  characterization. 

In  Chapter  4,  a  class  of  solution  techniques, 
generalized  product  form,  is  presented.  The  generalized 


product  form  approximation  yields  an  exact  solution  for  a 
queueing  network  satisfying  conditions  for  local  balance. 
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Most  of  the  existing  approximations  can  be  viewed  in  the 
framework  of  generalized  product  form-  Two  members  of  the 
generalized  product  form,  GPF  with  diffusion  and  GPF  with 
simulation,  are  presented.  These  approximations  use  the 
techniques  developed  in  Chapter  3  for  transition  process 
characterizations.  A  class  of  queueing  network  models  are 
considered  to  examine  the  accuracy  of  the  proposed 
approximations.  Both  of  the  proposed  approximations  show 
improved  accuracy  over  most  existing  approximations. 

In  Chapter  5,  the  accuracy  of  various  approximations 
for  a  family  of  parameters  is  examined.  In  the  absence  of 
any  analytic  expressions  for  error  bounds,  an  experimental 
approach  is  taken.  The  relationships  among  different 
parameters  and  the  errors  in  performance  measures,  e.g., 
mean  queue  lengths  and  utilizations,  are  explained. 


6.2  Directions  For  Future  Research 

The  work  presented  in  this  thesis  is  only  an  initial 
step  towards  understanding  and  improving  various,  seemingly 
unrelated,  approximations.  Further  efforts  are  needed  in 
this  direction  for  many  other  kinds  of  networks  and  many 
other  approximations.  In  this  section  we  briefly  outline 
some  of  the  directions  which  future  research  might  take. 

Diffusion  approximations  for  multiple  class  networks 
are  known  only  for  open  systems  [Gelenbe  and  Pujolle  1977]. 
To  extend  this  result  for  closed  networks  with  multiple 
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classes  the  generalized  product  form  framework  may  be 
applicable-  The  major  difference  between  solving  open  and 
closed  networks  is  the  availability  of  the  absolute  arrival 
rates  in  open  networks-  An  iterative  approach  for  multiple 
classes/  similar  to  GPF ,  can  be  used  to  obtain  the  absolute 
arrival  rates  in  closed  networks.  The  throughput  constraint 
in  this  case  has  to  be  satisfied  for  each  class 
individually. 

Our  experiments  for  determining  the  error  of  various 
approximations  have  been  limited  to  the  central  server, 
single  class  queueing  network  models  of  computer  systems. 
This  work  can  be  extended  in  two  directions:  general 
queueing  network  topologies  and  multiple  class  models.  We 
have  concentrated  on  the  central  server  model  since  it  has 
been  extensively  used  in  modeling  computer  systems  and  we 
were  primarily  interested  in  queueing  network  models  of 
computer  systems.  It  would,  however,  be  worthwhile  to 
analyze  the  approximations  for  general  topologies. 

Based  on  our  study  we  feel  that  the  results  obtained 
for  varying  the  coefficient  of  variation  of  the  service  time 
distribution  in  central  server  models  hold  for  non-central 
server  networks  as  well.  For  the  multiprogramming  level,  N, 
in  a  non-central  server  network,  the  error  behaviours  of 
various  approximations  should  be  similar  to  those  for 
the  central  server  model.  That  is,  if  N  is  greater  than  the 
system  saturation  level,  Nmax,  errors  in  utilizations  of  the 
servers  as  well  as  in  the  mean  queue  lengths  at  the  servers 


-T42- 


should  decrease.  For  smaller  values  of  N  the  error 
behaviour  needs  further  investigation. 

The  error  analysis  done  in  central  server  models  for 
various  approximations  when  the  number  of  servers  is  varied 
is  harder  to  generalize  for  non-central  server  networks.  In 
general,  the  effect  of  the  number  of  servers  on  the  accuracy 
of  an  approximation  will  depend  on  the  network  topology. 

With  an  extension  to  include  more  than  one  class  of 
customers  the  model  will  be  complex.  For  these  complex 
models  it  will  be  very  unlikely  that  exact  efficient 
solutions  will  be  available.  Hence,  for  error  analysis  of 


approximations 

for  the  multiple  class  case,  one  may  have  to 

consider  using 

simulation,  which  requires  significantly  more 

computer  time 

and  resources.  Further  efforts  in  multiple 

class  networks 

and  networks  with  priority  disciplines  are 

needed  to  understand  various  approximations. 

GPF  with  simulation  gives  very  accurate  results.  The 


time  required 

to  attain  a  given  level  of  accuracy  by 

simulating  the 

entire  network  is  significantly  greater  than 

that  required 

to  simulate  an  individual  centre,  times  the 

number  of  centres.  However,  GPF  with  simulation  requires 
substantially  more  computation  than  GPF  with  diffusion- 
selection  between  the  two  requires  a  trade-off  between 
accuracy  and  computational  expense. 

To  test  the  accuracy  of  a  new/old  approximation  there 
are  basically  two  methods:  obtain  an  analytical  error  bound 
or  test  the  accuracy  for  a  wide  range  of  network  topologies 
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and  parameters.  Both  of  these  methods  need  further 
investigation.  The  importance  of  the  first  technique  is 
well  recognized  and  efforts  are  being  made  to  obtain  error 
bounds  for  various  approximations  (although  with  very  little 
success).  We  feel  that  the  second  problem  needs  more 
attention.  Can  a  "critical"  set  of  queueing  networks  be 
constructed  to  test  the  accuracy  of  the  different 
approximations.  By  critical,  we  mean  a  set  that  includes 
topologies  and  parameters  for  queueing  networks  for  which 
most  approximations  have  their  worst  performance.  Creating 
such  a  set  should  be  very  helpful  for  the  error  analysis  of 
the  various  approximations. 
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APPENDIX  A 

A.l  Derivation  of  Equation  (3.1) 

Consider  the  distribution  of  the  number  of  events  routed  to 
stream  y2  between  each  pair  of  successive  events  routed  to  stream 
yl.  This  distribution  is: 

Prob  [k  events  of  stream  y2  between  two  events  to  stream  yl] 

=  P2k*Pl  k  =  0,1,2,.... 

k 

Therefore,  with  probability  P2  *P1,  an  interevent  time  in  stream  yl 
is  composed  of  the  sum  of  k  +  1  interevent  times  of  stream  x.  That 
is,  the  interevent  time  in  stream  yl  is  a  sample  from  the  (k+l)-fold 
convolution  of  the  distribution  of  stream  x  interevent  times. 
Expressions  for  the  mean  and  second  moment  of  the  distribution  of 
interevent  times  in  stream  yl  are  derived  below. 

Mean  (Myi ) :  ^ 

Myl  =  ^  Pl*P2k*E[sample  from  (K  +  l)-fold  convolution 
k=0  from  Fx(t)] 

00 

»  2  Pl*P2k*(-l)  [Lx(s)]k+1|s=0 

k=0 

where  L  (s)  is  the  Laplace  transform  of  F  (t). 
x  x 

Myl  =  £  -Pl*P2k  (k+1)  [Lx(s)]k|s=0*[L‘x(s)]|s=0 
k=0 


=  Pl*E[x]  2  (k+l)p2,< 
k=0 
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Mx 

PI 


Myl 


Mx 
"  PI 


Similarly,  My 2  =  . 


Second  Moment  (Syi): 


[using  the  fact  that  PI  =  1  -  P2] 


Syl  =  ^  Pl*P2k*E  [square  of  the  sample  from  (k+l)-fold 

convolution  from  Fx(t)] 


k=0 


£  Pl*P2k  [(Lx(s))k+1]| 

k=0  ds 


*  ^2  Pl*P2k*(k+l)[(Lx(s))kL"x(s)+(k)*(Lx(s))k'1* 


k=0 


U;(x)ni„ 


s=0 


Sx  ,  2*P2*(Mx)2 
'  P1  (PI)2 

CVyl  =  -  1 

( My 1 )  ^ 

=  Sx  +  2*P2*(Mx  )2  _  x 

PI* (My l)2  (Pl)2*(Myl)2 

=  Pl*(CVx+l)  +  2*P2-1 

=  Pl*CVx  +  PI  +  2  ( 1  -  PI )  -  1  CP2  =  1-Pl] 

=  Pl*CVx  +  1  -  PI 
=  1  +  Pl*(CVx-l) 


Similarly , 

CVy2  =  1  +  P2(CVx-l) . 


This  completes  the  derivation  of  equation  (3.1) 
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A.2  Proof  of  "the  first  n  moments  of  the  process  x  are  sufficient  to 

calculate  the  first  n  moments  of  process  yi . " 


i,L 

Let  the  rr  moments  of  interevent  times  in  processes  x  and  yi 

n  ,  n 
be  Mx  and  Myi . 


Myl  =  }  '  Pi*P2k*E  [n^1  power  of  sample  from  (k+l)-fold  convolution 
“  of  Fx ( t ) ] 


^Pl‘P2k%Lx(s))ktl]|s=0 
k=0  as 


CD 

=  y  ]  Pl*P2k[(k+]){(Lv(s)  )k*L^(s)  +  terms  involving  derivatives 

of  order  (n-1)  or  lower}] |s=0 

t 

where  L "  (s)  is  the  nth  derivative  of  L  (s). 

A  A 


n  4-  u 

Thus,  Myl  is  represented  as  the  sum  of  the  n  and  lower  moments 

of  Fx(t ) . 


This  completes  the  proof. 
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A.3  Splitting  into  More  Than  Two  Streams: 

Suppose  we  want  to  split  a  stream  into  n  substreams  with  probabil¬ 
ities  PI,  P2,  . .  ,  Pn  (PI  +  P2  +  ..  +  Pn  =  1).  The  n  substreams  are  divided 
into  two  groups  and  the  result  obtained  for  splitting  into  two  is  used 
to  calculate  the  means  and  the  CV's  of  the  two  substreams.  These  sub¬ 
streams  are  further  split  into  two  groups.  The  splitting  of  these 
substreams  continues  until  individual  substreams  are  obtained,  i.e.,  all 
the  n  resulting  streams  are  obtained.  There  are  many  possible  ways  by 
which  the  desired  splitting  can  be  obtained.  We  show  here  that  the  order 
in  which  these  groupings  are  done  has  no  effect  on  the  final  result. 

Let  us  assume  .that  we  want  to  calculate  the  mean  and  CV  for  the 
ith  (arbitrary  but  fixed)  resulting  stream  (Myi ,  CVyi) 

Let 

Ek  =  {Pj  |  Pj's  are  combined  with  the  designated  Pi  at  the  kth 
stage  of  splitting) 

Clearly,  Ek  3  Ek+1 

Since  there  are  finite  numbers  of  Pj's  ,  there  exists  an  m  such  that 
Em  =  (Pi  ,  Pm2) 

Pm2  may  be  any  of  the  PI,  P2,  ..,  Pi-1,  Pi+1,  ..,  Pn.  Let  CVk  and 
Mk  be  the  coefficient  of  variation  and  the  mean  of  the  inter- event 
time  distribution  corresponding  to  group  Ek  .  Also,  let  Sk  be  the 
sum  of  the  probabilities  in  group  Ek  .  By  equation  (3.1): 

Ml  =  Mx/Sl 

CV1  =  1  +  SI  *(CVx  -  1) 
or, 

(CV1-1)  =  SI  *(CVx  -  1)  0) 

The  splitting  of  El  into  two  groups  yields  E2  and  the  corresponding 


probability  is  S2/S1  . 
Therefore, 
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M2  =  Ml  /(S2/S1 )  =  Mx/S2 
(CV2-1 )  =  (S2/S1 ) (CV1 -1 ) 
Substituting  for  (CVl-l)from  (1),  we  get: 

(CV2-1 )  =  S2*(CVx-l ) 

By  induction  on  k  we  can  show  that 
Mk  =  Mx/Sk 
(CVk-1)  =  Sk*(CVx-l ) 

Therefore, 


Mm  =  Mx/Sm 

(CVm-1 )  =  Sm*(CVx-l ) 

Em  contains  only  two  elements,  thus  (3.1)  can  be  applied.  The 
probability  corresponding  to  the  ith  resulting  stream  is  Pi /Sm  . 
Therefore, 

Myi  =  Mm/ ( Pi /Sm) 

=  Mx/Pi 

and, 


( CVyi-1 )  =  ( Pi/Sm)( CVm-1 ) 


=  Pi*(CVx-l ) 


Therefore, 


Myi  =  Mx 

and  CVyi  =  1  +  Pi*(CVx-l) 

Since  arbitrary  grouping  was  chosen  and  this  does  not  affect  the 
final  result,  we  have  proved  that  all  grouping  would  yield  the  same 


result. 
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A.4  Merging  of  More  Than  Two  Streams: 

We  want  to  prove  that  merging  of  more  than  two  streams  yields  the  same 
result  independent  of  the  ordering  in  which  the  merging  of  two  streams  is 
used.  We  will  prove  this  by  induction.  Suppose  we  are  interested  in 
merging  n  streams  having  the  mean  inter-event  time  equal  to  Mxi  and 
the  C V  of  the  inter-event  event  time  distribution  equal  to  CVxi  .  The 
resulting  stream  has  mean  My  and  CV  equal  to  CVy  . 

For  n  =  2  there  is  only  one  way  to  merge  and  therefore  the  result 
is  true. 

For  n  =  3  there  are  three  different  ways  to  merge  in  groups  of  two: 


It  can  be  easily  verified  that  using  the  result  (3.3)  (merging  of  two 
processes)  all  three  ways  give  same  result.  Namely 
(1/My)  =  (1/Mxl)  +  (1/Mx2)  +  (1/Mx3) 

and  CVy  =  1  +  '£  (&V  (CVxi-l ) 

i=l  \  Mxi  J 


Now,  let  us  suppose  that  the  result  is  true  (i.e.,  the  ordering  in  which 
we  group  them  has  no  effect)  for  n  =  k  .  We  want  to  show  that  it  is 
true  for  n  =  k+1  . 

In  merging  k+1  processes  at  the  final  stage  we  have  the  following 
situation  : 

(Myl,CVyl)  1st  component^'  merged  streams) 


(My2}cVy2)  2nd  component (j  merged  streams) 


where  i  +  j  =  k  +  1 


and 


i  ,  j  >  1 


(I) 
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=*  i  and  j  both  are  less  than  equal  to  k  .  Thus,  no  matter  what  order¬ 
ing  is  used  to  derive  1st  and  2nd  component  by  the  induction  hypothesis 
the  results  are  the  same.  This  is  true  for  all  i  ,  j  combinations  such 
that  condition  (I)  is  satisfied. 

That  is, 

Myl  =  (1/Mxl)  +  (1/MX2)  +  ..  +  (1/Mxi) 

My2  =  (T/M  xD  +  (1/Mx2)  +  ..  +  (1/Mxj) 

CVyl  =  1  +  V  /  M_\  2  (CVxh-1 ) 

—  Mxh  / 

CVy2  =  1  +  ib/! )2  (CVxh-1) 
h=l  '  Mxh  / 

where  {(Mxh  ,  CVxh)}  and  {(Mxh  ,  CVxh)} 


are  reordering  of  {(Mxh  ,  CVxh)}. 


Therefore, 


(1/My)  =  (1/M.yl)  +  (1/My2) 

=  (1/Mxl)  +  (1/Mx2)  +  ..  +  (1/Mxi )  +  (1/Mxl)  +  ..  +  (1/Mxj) 
By  reordering  this  can  be  rewritten  as 

(1/My)  =  (1/Mxl)  +  ..  +  (1/Mxk+l) 


2 


2 


•,  2 


(CVxh-1) 


j 


(CVxh-1) 
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=  -1  +4-/My_\2  (cVxh-1 )  + 
hflVMxh  ) 


(CVxH-l 


By  rearranging  this  can  be  written  as 

CVy  =  1  +yV  \2  (CVxh-1) 

htiVMxh; 


This  proves  the  induction  hypothesis  and  thus  the  result. 
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APPENDIX  B 

B.l  Derivation  of  Equation  (3.5) 

Server  1  represents  the  CPU  and  servers  2,3..,M  represent  I/O 

"hh  * 

devices.  The  arrival  process  at  the  i  server  (i>2)  is  obtained 
by  splitting  the  departure  processes  from  the  first  server.  There¬ 
fore,  by  (3.1) 

CVai  =  1  +  Pi* (CVdl  -  1)  i  =  2,3,..,M  (B.l) 

Given  an  arrival  process  at  the  ith  server  (i  =  2,3,..,M),  the  C V 

J.L. 

of  the  departure  process  at  the  i&  server  is  given  by  (using  (3.4)): 
CVdl  =  x  +  (^-)  *(CVsi  -  1)  +  (l-({jffj  )*(CVa1  -  1)  (B.2) 

Substituting  for  CVai  from  (B.l),  (B.2)  yields 

/  2  2 

CVdi  =  1  +  (Hf)  *(CVsi  -  1)  +  (1  “  (rIt)  )*Pi*(CVdl  -  1)  (B.3) 

For  the  central  server,  server  1,  the  CV  of  the  arrival  and  departure 

processes  are  given  by  (using  (3.3)  and  (3.4)): 

M  2 

CVai  =  1  +  2  (H})  *(CVdi  -  1)  (B.4) 

i=2  '  7 

CVdl  =  1  +  (sr)2*(CVsl  -  1)  +  (1  -  (|jff)  )*<CVal  - 


(B.5) 
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Combining  (B.3)  and  (B.5)  gives: 

CVdl  -  1  +  feV*(cVsi  -  1)  +  (1 

d-ftl )2)(CVal  -  «] 


(CVsl  -  1)  + 


Substituting  for  CVdi  from  (B.6)  in  (B.4),  we  have: 

M  >  .  o  /  .  /..  ,\2 


Cval  =  1  + 


or. 


CVal  =  1  + 


i 


M 


)*Pi- 


Mal\ 

Mslj 


*(CVsl  -  1)  +  (1 


/MalV 


)*(CVal  -  1) 


E(kt)  m)  -  1)  ♦  (l-(|i)  )*Pi*(|l)  *(Cvsl  . 

i=2 


1) 


i=2 


However,  in  a  central  server  model, 

-  Mal 
Mai  -  ■w 


(B.6) 


(B.7) 


Therefore , 


This  completes  the  derivation  of  equation  (3.5) 
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B.2  Derivation  of  Equation  3.6 


iiL 

The  CV  of  the  departure  process-  at  the  itn  server  is  given  by 
(using  3.4)): 

.  .  2  2 

CVdi  =  1  +  ({j|}-)  *(CVs1  -  1)  +  (1  -  (j|ij  )*(CVa1  -  1)  (B.8) 


If  Pji  is  the  transition  probability  from  server  j  to  i  and  CVdji  is 
the  coefficient  of  variation  of  the  part  for  the  transition  process 
that  goes  from  server  j  to  i ,  then  by  (3.1) 

CVdji  =  1  +  Pji*(CVdj  -  1)  (B.9) 


But, 


CVai  =  1  + 


j=l 


*( CVdji  -  1) 


(B. 10) 


(SaJ*  Pj1*)  rePresents  the  Proportion  of  arrivals  from  server  j  to 
the  total  arrivals  at  server  i 


Therefore,  by  (B.9)  and  (B.10) 


CVai  =  1  + 


i 

j=l 


i  *(CVdj  -  1) 


(B. 11) 


Substituting  for  CVai  from  (B.ll)  in  (B.8)  yields 


2 

CVdi  =  1  +  (0)  *(CVsi  -  1)  +  (1 


9  M 

/ Mai Vn  Y' 

\Msi  VMaj/ 

•  i  '  ' 

J=1 


"*Pji3*(CVdj  -  1)] 


This  completes  the  derivation  of  equation  (3.6) 
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B.3  Derivation  of  Equation  (3.7) 


Equation  (13)  of  Gelenbe  and  Pujolle's  [1976]  paper  in  the 
notation  of  chapter  3  reduces  to: 

L\2 


CVdi 


1  -  (if)  *<CVs1  +  D  + 


^-(ifM2#) + 


1  +  Mai* 


ZftCVdi  -  l)*Pji  ♦  1)  *(p-)J] 


This  can  be  rewritten  as 

cvdi  =  1  +  (lf)2*(CVsi  '  D  +  U 


M 

'(Slf))*['1  +  Mai*S((CVdj  -  1)* 


or, 


CVdi  =  1  + 


*  »*fe)]*  [-!  * 2 (Si)  ‘"-(s)®#)*11 

*■  - 

2  M 

(If)  *<CVs1  - 1) +  a  -  (if^^E®*^'* 

j=i 


(CVdj  -  1) 


>] 


or, 


CVdi  =  1  +  (Sff)2*(CVsi  -  1)  +  (1 
M 


- 

\Msi  / 


[S(8r) 

Lj  =  l 


.2 


i*Pji  *(CVdj  -  1) 


-] 


This  completes  the  derivation  of  equation  (3.7) 
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appendix  c 


C.l  Outline  for  the  Proof  of  Lemma  5.1 


In  the  following,  we  briefly  outline  the  proof  of  lemma 
5.1.  For  a  more  formal  treatment  the  reader  is  referred  to 
Feller's  book  [1973].  For  waiting  time  Wi  to  the  nearest 
renewal  epoch  in  the  ith  process  we  have: 


P  [Wi  4:  t] 


(1/Mi) 


Fi  (y) )dy 


Since  each  individual  renewal  process  is  rare  (formally 
known  as  infinitessimal  renewal  process),  for  a  fixed  i, 
Fi  (y)  is  small  and  Mi  is  large. 

Therefore , 


P [Wi  ^  t]  =  (t/Mi) 


But , 

P  [W  ^  t]  =  P  [min  (Wi)  <  t] 

Therefore , 

P[W  >  t]  =  P  [all  Wi  >  t] 

= (l-(t/Ml) ) * (l-(t/M2) ) *. . .* (1- (t/Mn) ) 

As  n  approaches  infinity  {( 1/M  1 )  +  ( 1/M  2 )+....+( 1/Mn )  } 

approaches  1/M.  Therefore,  as  n  approaches  infinity 


P[W  >  t]  =  exp ( -t*M) 
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Thus,  the  superposition  of  rare  renewal  processes  approaches 
a  Poisson  stream. 


C.2  Outline  for  the  Proof  of  Lemma  5.2 

For  a  fixed  n,  we  have  shown  in  Chapter  3  (Equation 
(3.3)  )  that 

CVi  =  1  +  Pi(CV-l)  i=l , 2, . . . , n 

As  n  increases  all  Pi  approach  zero  uniformly.  Thus  the 
second  term  on  the  right  hand  side  of  the  above  equation  can 
be  made  arbitrarily  close  to  zero  for  sufficiently  large  n. 
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